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();
48 fTrack[0]=fTrack[1]=-1;
49 fGhostChi2Cut = 1e6; // nothing done by default
53 for(Int_t i=0; i<100; i++) {
54 for (Int_t j=0; j<2; j++) {
58 fRawClusters = new TClonesArray("AliMUONRawCluster",1000);
61 //____________________________________________________________________________
62 AliMUONClusterFinderVS::~AliMUONClusterFinderVS()
64 // Reset tracks information
67 fRawClusters->Delete();
72 AliMUONClusterFinderVS::AliMUONClusterFinderVS(const AliMUONClusterFinderVS & clusterFinder):TObject(clusterFinder)
74 // Protected copy constructor
76 AliFatal("Not implemented.");
78 //____________________________________________________________________________
79 void AliMUONClusterFinderVS::ResetRawClusters()
81 // Reset tracks information
83 if (fRawClusters) fRawClusters->Clear();
85 //____________________________________________________________________________
86 void AliMUONClusterFinderVS::Decluster(AliMUONRawCluster *cluster)
88 // Decluster by local maxima
89 SplitByLocalMaxima(cluster);
91 //____________________________________________________________________________
92 void AliMUONClusterFinderVS::SplitByLocalMaxima(AliMUONRawCluster *c)
94 // Split complex cluster by local maxima
97 fInput->SetCluster(c);
99 fMul[0]=c->GetMultiplicity(0);
100 fMul[1]=c->GetMultiplicity(1);
103 // dump digit information into arrays
108 for (cath=0; cath<2; cath++) {
111 for (i=0; i<fMul[cath]; i++) {
113 fDig[i][cath]=fInput->Digit(cath, c->GetIndex(i, cath));
115 fIx[i][cath]= fDig[i][cath]->PadX();
116 fIy[i][cath]= fDig[i][cath]->PadY();
118 fQ[i][cath] = fDig[i][cath]->Signal();
119 // pad centre coordinates
121 GetPadC(fInput->DetElemId(), fIx[i][cath], fIy[i][cath], fX[i][cath], fY[i][cath], fZ[i][cath]);
122 } // loop over cluster digits
124 } // loop over cathodes
130 // Initialise and perform mathieson fits
131 Float_t chi2, oldchi2;
132 // ++++++++++++++++++*************+++++++++++++++++++++
133 // (1) No more than one local maximum per cathode plane
134 // +++++++++++++++++++++++++++++++*************++++++++
135 if ((fNLocal[0]==1 && (fNLocal[1]==0 || fNLocal[1]==1)) ||
136 (fNLocal[0]==0 && fNLocal[1]==1)) {
137 // Perform combined single Mathieson fit
138 // Initial values for coordinates (x,y)
140 // One local maximum on cathodes 1 and 2 (X->cathode 2, Y->cathode 1)
141 if (fNLocal[0]==1 && fNLocal[1]==1) {
142 fXInit[0]=c->GetX(1);
143 fYInit[0]=c->GetY(0);
144 // One local maximum on cathode 1 (X,Y->cathode 1)
145 } else if (fNLocal[0]==1) {
146 fXInit[0]=c->GetX(0);
147 fYInit[0]=c->GetY(0);
148 // One local maximum on cathode 2 (X,Y->cathode 2)
150 fXInit[0]=c->GetX(1);
151 fYInit[0]=c->GetY(1);
153 AliDebug(1,"cas (1) CombiSingleMathiesonFit(c)");
154 chi2=CombiSingleMathiesonFit(c);
155 // Int_t ndf = fgNbins[0]+fgNbins[1]-2;
156 // Float_t prob = TMath::Prob(Double_t(chi2),ndf);
157 // prob1->Fill(prob);
158 // chi2_1->Fill(chi2);
160 AliDebug(1,Form(" chi2 %f ",chi2));
162 c->SetX(0, fXFit[0]);
163 c->SetY(0, fYFit[0]);
171 c->SetX(0, fSeg2[0]->GetAnod(fInput->DetElemId(), c->GetX(0)));
172 c->SetX(1, fSeg2[1]->GetAnod(fInput->DetElemId(), c->GetX(1)));
174 // c->SetDetElemId(fInput->DetElemId());
175 // If reasonable chi^2 add result to the list of rawclusters
178 // If not try combined double Mathieson Fit
180 AliDebug(1," MAUVAIS CHI2 !!!\n");
181 if (fNLocal[0]==1 && fNLocal[1]==1) {
182 fXInit[0]=fX[fIndLocal[0][1]][1];
183 fYInit[0]=fY[fIndLocal[0][0]][0];
184 fXInit[1]=fX[fIndLocal[0][1]][1];
185 fYInit[1]=fY[fIndLocal[0][0]][0];
186 } else if (fNLocal[0]==1) {
187 fXInit[0]=fX[fIndLocal[0][0]][0];
188 fYInit[0]=fY[fIndLocal[0][0]][0];
189 fXInit[1]=fX[fIndLocal[0][0]][0];
190 fYInit[1]=fY[fIndLocal[0][0]][0];
192 fXInit[0]=fX[fIndLocal[0][1]][1];
193 fYInit[0]=fY[fIndLocal[0][1]][1];
194 fXInit[1]=fX[fIndLocal[0][1]][1];
195 fYInit[1]=fY[fIndLocal[0][1]][1];
198 // Initial value for charge ratios
201 AliDebug(1,"\n cas (1) CombiDoubleMathiesonFit(c)\n");
202 chi2=CombiDoubleMathiesonFit(c);
203 // Int_t ndf = fgNbins[0]+fgNbins[1]-6;
204 // Float_t prob = TMath::Prob(chi2,ndf);
205 // prob2->Fill(prob);
206 // chi2_2->Fill(chi2);
208 // Was this any better ??
209 AliDebug(1,Form(" Old and new chi2 %f %f ", oldchi2, chi2));
210 if (fFitStat!=0 && chi2>0 && (2.*chi2 < oldchi2)) {
212 // Split cluster into two according to fit result
215 AliDebug(1,"Do not Split");
221 // +++++++++++++++++++++++++++++++++++++++
222 // (2) Two local maxima per cathode plane
223 // +++++++++++++++++++++++++++++++++++++++
224 } else if (fNLocal[0]==2 && fNLocal[1]==2) {
226 // Let's look for ghosts first
228 Float_t xm[4][2], ym[4][2];
229 Float_t dpx, dpy, dx, dy;
230 Int_t ixm[4][2], iym[4][2];
231 Int_t isec, im1, im2, ico;
233 // Form the 2x2 combinations
234 // 0-0, 0-1, 1-0, 1-1
236 for (im1=0; im1<2; im1++) {
237 for (im2=0; im2<2; im2++) {
238 xm[ico][0]=fX[fIndLocal[im1][0]][0];
239 ym[ico][0]=fY[fIndLocal[im1][0]][0];
240 xm[ico][1]=fX[fIndLocal[im2][1]][1];
241 ym[ico][1]=fY[fIndLocal[im2][1]][1];
243 ixm[ico][0]=fIx[fIndLocal[im1][0]][0];
244 iym[ico][0]=fIy[fIndLocal[im1][0]][0];
245 ixm[ico][1]=fIx[fIndLocal[im2][1]][1];
246 iym[ico][1]=fIy[fIndLocal[im2][1]][1];
250 // ico = 0 : first local maximum on cathodes 1 and 2
251 // ico = 1 : fisrt local maximum on cathode 1 and second on cathode 2
252 // ico = 2 : second local maximum on cathode 1 and first on cathode 1
253 // ico = 3 : second local maximum on cathodes 1 and 2
255 // Analyse the combinations and keep those that are possible !
256 // For each combination check consistency in x and y
259 Float_t dr[4] = {1.e4, 1.e4, 1.e4, 1.e4};
262 // In case of staggering maxima are displaced by exactly half the pad-size in y.
263 // We have to take into account the numerical precision in the consistency check;
266 for (ico=0; ico<4; ico++) {
267 accepted[ico]=kFALSE;
268 // cathode one: x-coordinate
269 isec=fSeg2[0]->Sector(fInput->DetElemId(), ixm[ico][0], iym[ico][0]);
270 dpx=fSeg2[0]->Dpx(fInput->DetElemId(), isec)/2.;
272 dx=TMath::Abs(xm[ico][0]-xm[ico][1]);
273 // cathode two: y-coordinate
275 isec=fSeg2[1]->Sector(fInput->DetElemId(), ixm[ico][1], iym[ico][1]);
276 dpy=fSeg2[1]->Dpy(fInput->DetElemId(), isec)/2.;
278 dy=TMath::Abs(ym[ico][0]-ym[ico][1]);
279 AliDebug(2,Form("\n %i %f %f %f %f %f %f \n", ico, ym[ico][0], ym[ico][1], dy, dpy, dx, dpx ));
280 if ((dx <= dpx) && (dy <= dpy+eps)) {
283 dr[ico] = TMath::Sqrt(dx*dx+dy*dy);
287 accepted[ico]=kFALSE;
290 AliDebug(1,Form("\n iacc= %d:\n", iacc));
292 if (accepted[0] && accepted[1]) {
293 if (dr[0] >= dr[1]) {
300 if (accepted[2] && accepted[3]) {
301 if (dr[2] >= dr[3]) {
308 // eliminate one candidate
312 for (ico=0; ico<4; ico++) {
313 if (accepted[ico] && dr[ico] > drmax) {
319 accepted[icobad] = kFALSE;
325 AliDebug(1,Form("\n iacc= %d:\n", iacc));
327 AliDebug(1,"\n iacc=2: No problem ! \n");
328 } else if (iacc==4) {
329 AliDebug(1,"\n iacc=4: Ok, but ghost problem !!! \n");
330 } else if (iacc==0) {
331 AliDebug(1,"\n iacc=0: I don't know what to do with this !!!!!!!!! \n");
334 // Initial value for charge ratios
335 fQrInit[0]=Float_t(fQ[fIndLocal[0][0]][0])/
336 Float_t(fQ[fIndLocal[0][0]][0]+fQ[fIndLocal[1][0]][0]);
337 fQrInit[1]=Float_t(fQ[fIndLocal[0][1]][1])/
338 Float_t(fQ[fIndLocal[0][1]][1]+fQ[fIndLocal[1][1]][1]);
340 // ******* iacc = 0 *******
341 // No combinations found between the 2 cathodes
342 // We keep the center of gravity of the cluster
347 // ******* iacc = 1 *******
348 // Only one combination found between the 2 cathodes
350 // Initial values for the 2 maxima (x,y)
352 // 1 maximum is initialised with the maximum of the combination found (X->cathode 2, Y->cathode 1)
353 // 1 maximum is initialised with the other maximum of the first cathode
360 } else if (accepted[1]){
366 } else if (accepted[2]){
372 } else if (accepted[3]){
379 AliDebug(1,"cas (2) CombiDoubleMathiesonFit(c)");
380 chi2=CombiDoubleMathiesonFit(c);
381 // Int_t ndf = fgNbins[0]+fgNbins[1]-6;
382 // Float_t prob = TMath::Prob(chi2,ndf);
383 // prob2->Fill(prob);
384 // chi2_2->Fill(chi2);
385 AliDebug(1,Form(" chi2 %f\n",chi2));
387 // If reasonable chi^2 add result to the list of rawclusters
392 // 1 maximum is initialised with the maximum of the combination found (X->cathode 2, Y->cathode 1)
393 // 1 maximum is initialised with the other maximum of the second cathode
400 } else if (accepted[1]){
406 } else if (accepted[2]){
412 } else if (accepted[3]){
419 AliDebug(1,"\n cas (2) CombiDoubleMathiesonFit(c)\n");
420 chi2=CombiDoubleMathiesonFit(c);
421 // Int_t ndf = fgNbins[0]+fgNbins[1]-6;
422 // Float_t prob = TMath::Prob(chi2,ndf);
423 // prob2->Fill(prob);
424 // chi2_2->Fill(chi2);
425 AliDebug(1,Form(" chi2 %f\n",chi2));
427 // If reasonable chi^2 add result to the list of rawclusters
431 //We keep only the combination found (X->cathode 2, Y->cathode 1)
432 for (Int_t ico=0; ico<2; ico++) {
434 AliMUONRawCluster cnew;
436 for (cath=0; cath<2; cath++) {
437 cnew.SetX(cath, Float_t(xm[ico][1]));
438 cnew.SetY(cath, Float_t(ym[ico][0]));
439 cnew.SetZ(cath, fZPlane);
440 cnew.SetMultiplicity(cath,c->GetMultiplicity(cath));
441 for (i=0; i<fMul[cath]; i++) {
442 cnew.SetIndex(i, cath, c->GetIndex(i,cath));
443 fSeg2[cath]->SetPad(fInput->DetElemId(), fIx[i][cath], fIy[i][cath]);
445 AliDebug(1,Form("\nRawCluster %d cath %d\n",ico,cath));
446 AliDebug(1,Form("mult_av %d\n",c->GetMultiplicity(cath)));
447 FillCluster(&cnew,cath);
449 cnew.SetClusterType(cnew.PhysicsContribution());
458 // ******* iacc = 2 *******
459 // Two combinations found between the 2 cathodes
461 // Was the same maximum taken twice
462 if ((accepted[0]&&accepted[1]) || (accepted[2]&&accepted[3])) {
463 AliDebug(1,"\n Maximum taken twice !!!\n");
465 // Have a try !! with that
466 if (accepted[0]&&accepted[3]) {
477 AliDebug(1,"\n cas (2) CombiDoubleMathiesonFit(c)\n");
478 chi2=CombiDoubleMathiesonFit(c);
479 // Int_t ndf = fgNbins[0]+fgNbins[1]-6;
480 // Float_t prob = TMath::Prob(chi2,ndf);
481 // prob2->Fill(prob);
482 // chi2_2->Fill(chi2);
486 // No ghosts ! No Problems ! - Perform one fit only !
487 if (accepted[0]&&accepted[3]) {
498 AliDebug(1,"\n cas (2) CombiDoubleMathiesonFit(c)\n");
499 chi2=CombiDoubleMathiesonFit(c);
500 // Int_t ndf = fgNbins[0]+fgNbins[1]-6;
501 // Float_t prob = TMath::Prob(chi2,ndf);
502 // prob2->Fill(prob);
503 // chi2_2->Fill(chi2);
504 AliDebug(1,Form(" chi2 %f\n",chi2));
508 // ******* iacc = 4 *******
509 // Four combinations found between the 2 cathodes
511 } else if (iacc==4) {
512 // Perform fits for the two possibilities !!
513 // Accept if charges are compatible on both cathodes
514 // If none are compatible, keep everything
519 AliDebug(1,"\n cas (2) CombiDoubleMathiesonFit(c)\n");
520 chi2=CombiDoubleMathiesonFit(c);
521 // Int_t ndf = fgNbins[0]+fgNbins[1]-6;
522 // Float_t prob = TMath::Prob(chi2,ndf);
523 // prob2->Fill(prob);
524 // chi2_2->Fill(chi2);
525 AliDebug(1,Form(" chi2 %f\n",chi2));
526 // store results of fit and postpone decision
527 Double_t sXFit[2],sYFit[2],sQrFit[2];
529 for (Int_t i=0;i<2;i++) {
539 AliDebug(1,"\n cas (2) CombiDoubleMathiesonFit(c)\n");
540 chi2=CombiDoubleMathiesonFit(c);
541 // ndf = fgNbins[0]+fgNbins[1]-6;
542 // prob = TMath::Prob(chi2,ndf);
543 // prob2->Fill(prob);
544 // chi2_2->Fill(chi2);
545 AliDebug(1,Form(" chi2 %f\n",chi2));
546 // We have all informations to perform the decision
547 // Compute the chi2 for the 2 possibilities
548 Float_t chi2fi,chi2si,chi2f,chi2s;
550 chi2f = (TMath::Log(fInput->TotalCharge(0)*fQrFit[0]
551 / (fInput->TotalCharge(1)*fQrFit[1]) )
552 / fInput->ChargeCorrel() );
554 chi2fi = (TMath::Log(fInput->TotalCharge(0)*(1-fQrFit[0])
555 / (fInput->TotalCharge(1)*(1-fQrFit[1])) )
556 / fInput->ChargeCorrel() );
557 chi2f += chi2fi*chi2fi;
559 chi2s = (TMath::Log(fInput->TotalCharge(0)*sQrFit[0]
560 / (fInput->TotalCharge(1)*sQrFit[1]) )
561 / fInput->ChargeCorrel() );
563 chi2si = (TMath::Log(fInput->TotalCharge(0)*(1-sQrFit[0])
564 / (fInput->TotalCharge(1)*(1-sQrFit[1])) )
565 / fInput->ChargeCorrel() );
566 chi2s += chi2si*chi2si;
568 // usefull to store the charge matching chi2 in the cluster
569 // fChi2[0]=sChi2[1]=chi2f;
570 // fChi2[1]=sChi2[0]=chi2s;
572 if (chi2f<=fGhostChi2Cut && chi2s<=fGhostChi2Cut)
574 if (chi2f>fGhostChi2Cut && chi2s>fGhostChi2Cut) {
580 if (chi2f<=fGhostChi2Cut)
582 if (chi2s<=fGhostChi2Cut) {
583 // retreive saved values
584 for (Int_t i=0;i<2;i++) {
595 } else if (fNLocal[0]==2 && fNLocal[1]==1) {
596 // ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
597 // (3) Two local maxima on cathode 1 and one maximum on cathode 2
598 // ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
600 Float_t xm[4][2], ym[4][2];
601 Float_t dpx, dpy, dx, dy;
602 Int_t ixm[4][2], iym[4][2];
603 Int_t isec, im1, ico;
605 // Form the 2x2 combinations
606 // 0-0, 0-1, 1-0, 1-1
608 for (im1=0; im1<2; im1++) {
609 xm[ico][0]=fX[fIndLocal[im1][0]][0];
610 ym[ico][0]=fY[fIndLocal[im1][0]][0];
611 xm[ico][1]=fX[fIndLocal[0][1]][1];
612 ym[ico][1]=fY[fIndLocal[0][1]][1];
614 ixm[ico][0]=fIx[fIndLocal[im1][0]][0];
615 iym[ico][0]=fIy[fIndLocal[im1][0]][0];
616 ixm[ico][1]=fIx[fIndLocal[0][1]][1];
617 iym[ico][1]=fIy[fIndLocal[0][1]][1];
620 // ico = 0 : first local maximum on cathodes 1 and 2
621 // ico = 1 : second local maximum on cathode 1 and first on cathode 2
623 // Analyse the combinations and keep those that are possible !
624 // For each combination check consistency in x and y
628 // In case of staggering maxima are displaced by exactly half the pad-size in y.
629 // We have to take into account the numerical precision in the consistency check;
633 for (ico=0; ico<2; ico++) {
634 isec=fSeg2[0]->Sector(fInput->DetElemId(), ixm[ico][0], iym[ico][0]);
635 dpx=fSeg2[0]->Dpx(fInput->DetElemId(), isec)/2.;
637 dx=TMath::Abs(xm[ico][0]-xm[ico][1]);
638 isec=fSeg2[1]->Sector(fInput->DetElemId(), ixm[ico][1], iym[ico][1]);
639 dpy=fSeg2[1]->Dpy(fInput->DetElemId(), isec)/2.;
641 dy=TMath::Abs(ym[ico][0]-ym[ico][1]);
642 AliDebug(2,Form("\n %i %f %f %f %f \n", ico, ym[ico][0], ym[ico][1], dy, dpy ));
643 if ((dx <= dpx) && (dy <= dpy+eps)) {
649 accepted[ico]=kFALSE;
657 // Initial value for charge ratios
658 fQrInit[0]=Float_t(fQ[fIndLocal[0][0]][0])/
659 Float_t(fQ[fIndLocal[0][0]][0]+fQ[fIndLocal[1][0]][0]);
660 fQrInit[1]=fQrInit[0];
662 if (accepted[0] && accepted[1]) {
664 fXInit[0]=0.5*(xm[0][1]+xm[0][0]);
666 fXInit[1]=0.5*(xm[0][1]+xm[1][0]);
670 chi23=CombiDoubleMathiesonFit(c);
679 } else if (accepted[0]) {
684 chi21=CombiDoubleMathiesonFit(c);
685 // Int_t ndf = fgNbins[0]+fgNbins[1]-6;
686 // Float_t prob = TMath::Prob(chi2,ndf);
687 // prob2->Fill(prob);
688 // chi2_2->Fill(chi21);
689 AliDebug(1,Form(" chi2 %f\n",chi21));
690 if (chi21<10) Split(c);
691 } else if (accepted[1]) {
696 chi22=CombiDoubleMathiesonFit(c);
697 // Int_t ndf = fgNbins[0]+fgNbins[1]-6;
698 // Float_t prob = TMath::Prob(chi2,ndf);
699 // prob2->Fill(prob);
700 // chi2_2->Fill(chi22);
701 AliDebug(1,Form(" chi2 %f\n",chi22));
702 if (chi22<10) Split(c);
705 if (chi21 > 10 && chi22 > 10 && chi23 > 10) {
706 // We keep only the combination found (X->cathode 2, Y->cathode 1)
707 for (Int_t ico=0; ico<2; ico++) {
709 AliMUONRawCluster cnew;
711 for (cath=0; cath<2; cath++) {
712 cnew.SetX(cath, Float_t(xm[ico][1]));
713 cnew.SetY(cath, Float_t(ym[ico][0]));
714 cnew.SetZ(cath, fZPlane);
715 cnew.SetMultiplicity(cath, c->GetMultiplicity(cath));
716 for (i=0; i<fMul[cath]; i++) {
717 cnew.SetIndex(i, cath, c->GetIndex(i, cath));
718 fSeg2[cath]->SetPad(fInput->DetElemId(), fIx[i][cath], fIy[i][cath]);
721 AliDebug(1,Form("\nRawCluster %d cath %d\n",ico,cath));
722 AliDebug(1,Form("mult_av %d\n",c->GetMultiplicity(cath)));
724 FillCluster(&cnew,cath);
726 cnew.SetClusterType(cnew.PhysicsContribution());
733 // +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
734 // (3') One local maximum on cathode 1 and two maxima on cathode 2
735 // +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
736 } else if (fNLocal[0]==1 && fNLocal[1]==2) {
737 Float_t xm[4][2], ym[4][2];
738 Float_t dpx, dpy, dx, dy;
739 Int_t ixm[4][2], iym[4][2];
740 Int_t isec, im1, ico;
742 // Form the 2x2 combinations
743 // 0-0, 0-1, 1-0, 1-1
745 for (im1=0; im1<2; im1++) {
746 xm[ico][0]=fX[fIndLocal[0][0]][0];
747 ym[ico][0]=fY[fIndLocal[0][0]][0];
748 xm[ico][1]=fX[fIndLocal[im1][1]][1];
749 ym[ico][1]=fY[fIndLocal[im1][1]][1];
751 ixm[ico][0]=fIx[fIndLocal[0][0]][0];
752 iym[ico][0]=fIy[fIndLocal[0][0]][0];
753 ixm[ico][1]=fIx[fIndLocal[im1][1]][1];
754 iym[ico][1]=fIy[fIndLocal[im1][1]][1];
757 // ico = 0 : first local maximum on cathodes 1 and 2
758 // ico = 1 : first local maximum on cathode 1 and second on cathode 2
760 // Analyse the combinations and keep those that are possible !
761 // For each combination check consistency in x and y
765 // In case of staggering maxima are displaced by exactly half the pad-size in y.
766 // We have to take into account the numerical precision in the consistency check;
770 for (ico=0; ico<2; ico++) {
771 accepted[ico]=kFALSE;
772 isec=fSeg2[0]->Sector(fInput->DetElemId(), ixm[ico][0], iym[ico][0]);
773 dpx=fSeg2[0]->Dpx(fInput->DetElemId(), isec)/2.;
775 dx=TMath::Abs(xm[ico][0]-xm[ico][1]);
776 isec=fSeg2[1]->Sector(fInput->DetElemId(), ixm[ico][1], iym[ico][1]);
777 dpy=fSeg2[1]->Dpy(fInput->DetElemId(), isec)/2.;
779 dy=TMath::Abs(ym[ico][0]-ym[ico][1]);
780 AliDebug(1,Form("\n %i %f %f %f %f \n", ico, ym[ico][0], ym[ico][1], dy, dpy ));
781 if ((dx <= dpx) && (dy <= dpy+eps)) {
784 AliDebug(1,Form("ico %d\n",ico));
788 accepted[ico]=kFALSE;
796 fQrInit[1]=Float_t(fQ[fIndLocal[0][1]][1])/
797 Float_t(fQ[fIndLocal[0][1]][1]+fQ[fIndLocal[1][1]][1]);
799 fQrInit[0]=fQrInit[1];
802 if (accepted[0] && accepted[1]) {
804 fYInit[0]=0.5*(ym[0][0]+ym[0][1]);
806 fYInit[1]=0.5*(ym[0][0]+ym[1][1]);
809 chi23=CombiDoubleMathiesonFit(c);
818 } else if (accepted[0]) {
823 chi21=CombiDoubleMathiesonFit(c);
824 // Int_t ndf = fgNbins[0]+fgNbins[1]-6;
825 // Float_t prob = TMath::Prob(chi2,ndf);
826 // prob2->Fill(prob);
827 // chi2_2->Fill(chi21);
828 AliDebug(1,Form(" chi2 %f\n",chi21));
829 if (chi21<10) Split(c);
830 } else if (accepted[1]) {
835 chi22=CombiDoubleMathiesonFit(c);
836 // Int_t ndf = fgNbins[0]+fgNbins[1]-6;
837 // Float_t prob = TMath::Prob(chi2,ndf);
838 // prob2->Fill(prob);
839 // chi2_2->Fill(chi22);
840 AliDebug(1,Form(" chi2 %f\n",chi22));
841 if (chi22<10) Split(c);
844 if (chi21 > 10 && chi22 > 10 && chi23 > 10) {
845 //We keep only the combination found (X->cathode 2, Y->cathode 1)
846 for (Int_t ico=0; ico<2; ico++) {
848 AliMUONRawCluster cnew;
850 for (cath=0; cath<2; cath++) {
851 cnew.SetX(cath, Float_t(xm[ico][1]));
852 cnew.SetY(cath, Float_t(ym[ico][0]));
853 cnew.SetZ(cath, fZPlane);
854 cnew.SetMultiplicity(cath, c->GetMultiplicity(cath));
855 for (i=0; i<fMul[cath]; i++) {
856 cnew.SetIndex(i, cath, c->GetIndex(i, cath));
857 fSeg2[cath]->SetPad(fInput->DetElemId(), fIx[i][cath], fIy[i][cath]);
859 AliDebug(1,Form("\nRawCluster %d cath %d\n",ico,cath));
860 AliDebug(1,Form("mult_av %d\n",c->GetMultiplicity(cath)));
861 FillCluster(&cnew,cath);
863 cnew.SetClusterType(cnew.PhysicsContribution());
870 // ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
871 // (4) At least three local maxima on cathode 1 or on cathode 2
872 // ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
873 } else if (fNLocal[0]>2 || fNLocal[1]>2) {
874 Int_t param = fNLocal[0]*fNLocal[1];
877 Float_t ** xm = new Float_t * [param];
878 for (ii=0; ii<param; ii++) xm[ii]=new Float_t [2];
879 Float_t ** ym = new Float_t * [param];
880 for (ii=0; ii<param; ii++) ym[ii]=new Float_t [2];
881 Int_t ** ixm = new Int_t * [param];
882 for (ii=0; ii<param; ii++) ixm[ii]=new Int_t [2];
883 Int_t ** iym = new Int_t * [param];
884 for (ii=0; ii<param; ii++) iym[ii]=new Int_t [2];
887 Float_t dpx, dpy, dx, dy;
890 for (Int_t im1=0; im1<fNLocal[0]; im1++) {
891 for (Int_t im2=0; im2<fNLocal[1]; im2++) {
892 xm[ico][0]=fX[fIndLocal[im1][0]][0];
893 ym[ico][0]=fY[fIndLocal[im1][0]][0];
894 xm[ico][1]=fX[fIndLocal[im2][1]][1];
895 ym[ico][1]=fY[fIndLocal[im2][1]][1];
897 ixm[ico][0]=fIx[fIndLocal[im1][0]][0];
898 iym[ico][0]=fIy[fIndLocal[im1][0]][0];
899 ixm[ico][1]=fIx[fIndLocal[im2][1]][1];
900 iym[ico][1]=fIy[fIndLocal[im2][1]][1];
906 AliDebug(1,Form("nIco %d\n",nIco));
907 for (ico=0; ico<nIco; ico++) {
908 AliDebug(1,Form("ico = %d\n",ico));
909 isec=fSeg2[0]->Sector(fInput->DetElemId(), ixm[ico][0], iym[ico][0]);
910 dpx=fSeg2[0]->Dpx(fInput->DetElemId(), isec)/2.;
912 dx=TMath::Abs(xm[ico][0]-xm[ico][1]);
913 isec=fSeg2[1]->Sector(fInput->DetElemId(), ixm[ico][1], iym[ico][1]);
914 dpy=fSeg2[1]->Dpy(fInput->DetElemId(), isec)/2.;
916 dy=TMath::Abs(ym[ico][0]-ym[ico][1]);
917 AliDebug(1,Form("dx %f dpx %f dy %f dpy %f\n",dx,dpx,dy,dpy));
918 AliDebug(1,Form(" X %f Y %f\n",xm[ico][1],ym[ico][0]));
919 if ((dx <= dpx) && (dy <= dpy)) {
922 AliMUONRawCluster cnew;
923 for (cath=0; cath<2; cath++) {
924 cnew.SetX(cath, Float_t(xm[ico][1]));
925 cnew.SetY(cath, Float_t(ym[ico][0]));
926 cnew.SetZ(cath, fZPlane);
927 cnew.SetMultiplicity(cath, c->GetMultiplicity(cath));
928 for (i=0; i<fMul[cath]; i++) {
929 cnew.SetIndex(i, cath, c->GetIndex(i, cath));
930 fSeg2[cath]->SetPad(fInput->DetElemId(), fIx[i][cath], fIy[i][cath]);
932 FillCluster(&cnew,cath);
934 cnew.SetClusterType(cnew.PhysicsContribution());
935 // cnew.SetDetElemId(fInput->DetElemId());
947 void AliMUONClusterFinderVS::FindLocalMaxima(AliMUONRawCluster* /*c*/)
949 // Find all local maxima of a cluster
950 AliDebug(1,"\n Find Local maxima !");
954 Int_t cath, cath1; // loops over cathodes
955 Int_t i; // loops over digits
956 Int_t j; // loops over cathodes
960 // counters for number of local maxima
961 fNLocal[0]=fNLocal[1]=0;
962 // flags digits as local maximum
963 Bool_t isLocal[100][2];
964 for (i=0; i<100;i++) {
965 isLocal[i][0]=isLocal[i][1]=kFALSE;
967 // number of next neighbours and arrays to store them
970 // loop over cathodes
971 for (cath=0; cath<2; cath++) {
972 // loop over cluster digits
973 for (i=0; i<fMul[cath]; i++) {
974 // get neighbours for that digit and assume that it is local maximum
978 fSeg2[cath]->Neighbours(fInput->DetElemId(), fIx[i][cath], fIy[i][cath], &nn, x, y);
980 isLocal[i][cath]=kTRUE;
981 isec = fSeg2[cath]->Sector(fInput->DetElemId(), fIx[i][cath], fIy[i][cath]);
982 a0 = fSeg2[cath]->Dpx(fInput->DetElemId(), isec)*fSeg2[cath]->Dpy(fInput->DetElemId(), isec);
984 // loop over next neighbours, if at least one neighbour has higher charger assumption
985 // digit is not local maximum
986 for (j=0; j<nn; j++) {
987 if (fHitMap[cath]->TestHit(x[j], y[j])==kEmpty) continue;
988 digt=(AliMUONDigit*) fHitMap[cath]->GetHit(x[j], y[j]);
990 isec=fSeg2[cath]->Sector(fInput->DetElemId(), x[j], y[j]);
991 a1 = fSeg2[cath]->Dpx(fInput->DetElemId(),isec)*fSeg2[cath]->Dpy(fInput->DetElemId(), isec);
993 if (digt->Signal()/a1 > fQ[i][cath]/a0) {
994 isLocal[i][cath]=kFALSE;
997 // handle special case of neighbouring pads with equal signal
998 } else if (digt->Signal() == fQ[i][cath]) {
999 if (fNLocal[cath]>0) {
1000 for (Int_t k=0; k<fNLocal[cath]; k++) {
1001 if (x[j]==fIx[fIndLocal[k][cath]][cath]
1002 && y[j]==fIy[fIndLocal[k][cath]][cath])
1004 isLocal[i][cath]=kFALSE;
1006 } // loop over local maxima
1007 } // are there already local maxima
1009 } // loop over next neighbours
1010 if (isLocal[i][cath]) {
1011 fIndLocal[fNLocal[cath]][cath]=i;
1014 } // loop over all digits
1015 } // loop over cathodes
1017 AliDebug(1,Form("\n Found %d %d %d %d local Maxima\n",
1018 fNLocal[0], fNLocal[1], fMul[0], fMul[1]));
1019 AliDebug(1,Form("\n Cathode 1 local Maxima %d Multiplicite %d\n",fNLocal[0], fMul[0]));
1020 AliDebug(1,Form(" Cathode 2 local Maxima %d Multiplicite %d\n",fNLocal[1], fMul[1]));
1025 if (fNLocal[1]==2 && (fNLocal[0]==1 || fNLocal[0]==0)) {
1026 Int_t iback=fNLocal[0];
1028 // Two local maxima on cathode 2 and one maximum on cathode 1
1029 // Look for local maxima considering up and down neighbours on the 1st cathode only
1031 // Loop over cluster digits
1035 for (i=0; i<fMul[cath]; i++) {
1036 isec=fSeg2[cath]->Sector(fInput->DetElemId(), fIx[i][cath],fIy[i][cath]);
1037 dpy=fSeg2[cath]->Dpy(fInput->DetElemId(), isec);
1038 dpx=fSeg2[cath]->Dpx(fInput->DetElemId(), isec);
1040 if (isLocal[i][cath]) continue;
1041 // Pad position should be consistent with position of local maxima on the opposite cathode
1042 if ((TMath::Abs(fX[i][cath]-fX[fIndLocal[0][cath1]][cath1]) > dpx/2.) &&
1043 (TMath::Abs(fX[i][cath]-fX[fIndLocal[1][cath1]][cath1]) > dpx/2.))
1046 // get neighbours for that digit and assume that it is local maximum
1047 isLocal[i][cath]=kTRUE;
1048 // compare signal to that on the two neighbours on the left and on the right
1049 // iNN counts the number of neighbours with signal, it should be 1 or 2
1053 for (fSeg2[cath]->FirstPad(fInput->DetElemId(), fX[i][cath], fY[i][cath], fZPlane, 0., dpy);
1054 fSeg2[cath]->MorePads(fInput->DetElemId());
1055 fSeg2[cath]->NextPad(fInput->DetElemId()))
1057 ix = fSeg2[cath]->Ix();
1058 iy = fSeg2[cath]->Iy();
1059 // skip the current pad
1060 if (iy == fIy[i][cath]) continue;
1062 if (fHitMap[cath]->TestHit(ix, iy)!=kEmpty) {
1064 digt=(AliMUONDigit*) fHitMap[cath]->GetHit(ix,iy);
1065 if (digt->Signal() > fQ[i][cath]) isLocal[i][cath]=kFALSE;
1067 } // Loop over pad neighbours in y
1069 if (isLocal[i][cath] && iNN>0) {
1070 fIndLocal[fNLocal[cath]][cath]=i;
1073 } // loop over all digits
1074 // if one additional maximum has been found we are happy
1075 // if more maxima have been found restore the previous situation
1076 AliDebug(1,Form("\n New search gives %d local maxima for cathode 1 \n",
1078 AliDebug(1,Form(" %d local maxima for cathode 2 \n",
1080 if (fNLocal[cath]>2) {
1081 fNLocal[cath]=iback;
1084 } // 1,2 local maxima
1086 if (fNLocal[0]==2 && (fNLocal[1]==1 || fNLocal[1]==0)) {
1087 Int_t iback=fNLocal[1];
1089 // Two local maxima on cathode 1 and one maximum on cathode 2
1090 // Look for local maxima considering left and right neighbours on the 2nd cathode only
1093 Float_t eps = 1.e-5;
1096 // Loop over cluster digits
1097 for (i=0; i<fMul[cath]; i++) {
1098 isec=fSeg2[cath]->Sector(fInput->DetElemId(), fIx[i][cath],fIy[i][cath]);
1099 dpx=fSeg2[cath]->Dpx(fInput->DetElemId(), isec);
1100 dpy=fSeg2[cath]->Dpy(fInput->DetElemId(), isec);
1103 if (isLocal[i][cath]) continue;
1104 // Pad position should be consistent with position of local maxima on the opposite cathode
1105 if ((TMath::Abs(fY[i][cath]-fY[fIndLocal[0][cath1]][cath1]) > dpy/2.+eps) &&
1106 (TMath::Abs(fY[i][cath]-fY[fIndLocal[1][cath1]][cath1]) > dpy/2.+eps))
1110 // get neighbours for that digit and assume that it is local maximum
1111 isLocal[i][cath]=kTRUE;
1112 // compare signal to that on the two neighbours on the left and on the right
1114 // iNN counts the number of neighbours with signal, it should be 1 or 2
1116 for (fSeg2[cath]->FirstPad(fInput->DetElemId(), fX[i][cath], fY[i][cath], fZPlane, dpx, 0.);
1117 fSeg2[cath]->MorePads(fInput->DetElemId());
1118 fSeg2[cath]->NextPad(fInput->DetElemId()))
1121 ix = fSeg2[cath]->Ix();
1122 iy = fSeg2[cath]->Iy();
1124 // skip the current pad
1125 if (ix == fIx[i][cath]) continue;
1127 if (fHitMap[cath]->TestHit(ix, iy)!=kEmpty) {
1129 digt=(AliMUONDigit*) fHitMap[cath]->GetHit(ix,iy);
1130 if (digt->Signal() > fQ[i][cath]) isLocal[i][cath]=kFALSE;
1132 } // Loop over pad neighbours in x
1134 if (isLocal[i][cath] && iNN>0) {
1135 fIndLocal[fNLocal[cath]][cath]=i;
1138 } // loop over all digits
1139 // if one additional maximum has been found we are happy
1140 // if more maxima have been found restore the previous situation
1141 AliDebug(1,Form("\n New search gives %d local maxima for cathode 1 \n",fNLocal[0]));
1142 AliDebug(1,Form("\n %d local maxima for cathode 2 \n",fNLocal[1]));
1143 AliDebug(1,Form("\n New search gives %d %d \n",fNLocal[0],fNLocal[1]));
1144 if (fNLocal[cath]>2) {
1145 fNLocal[cath]=iback;
1147 } // 2,1 local maxima
1151 void AliMUONClusterFinderVS::FillCluster(AliMUONRawCluster* c, Int_t flag, Int_t cath)
1154 // Completes cluster information starting from list of digits
1161 c->SetPeakSignal(cath,c->GetPeakSignal(0));
1163 c->SetPeakSignal(cath,0);
1170 c->SetCharge(cath,0);
1173 AliDebug(1,Form("\n fPeakSignal %d\n",c->GetPeakSignal(cath)));
1174 for (Int_t i=0; i<c->GetMultiplicity(cath); i++)
1176 dig= fInput->Digit(cath,c->GetIndex(i,cath));
1177 ix=dig->PadX()+c->GetOffset(i,cath);
1179 Int_t q=dig->Signal();
1180 if (!flag) q=Int_t(q*c->GetContrib(i,cath));
1181 // fprintf(stderr,"q %d c->fPeakSignal[ %d ] %d\n",q,cath,c->fPeakSignal[cath]);
1182 if (dig->Physics() >= dig->Signal()) {
1184 } else if (dig->Physics() == 0) {
1186 } else c->SetPhysics(i,1);
1189 AliDebug(2,Form("q %d c->fPeakSignal[cath] %d\n",q,c->GetPeakSignal(cath)));
1190 // peak signal and track list
1191 if (q>c->GetPeakSignal(cath)) {
1192 c->SetPeakSignal(cath, q);
1193 c->SetTrack(0,dig->Hit());
1194 c->SetTrack(1,dig->Track(0));
1195 c->SetTrack(2,dig->Track(1));
1196 // fprintf(stderr," c->fTracks[0] %d c->fTracks[1] %d\n",dig->fHit,dig->fTracks[0]);
1200 fSeg2[cath]->GetPadC(fInput->DetElemId(), ix, iy, x, y, z);
1204 c->AddCharge(cath, q);
1206 } // loop over digits
1207 AliDebug(1," fin du cluster c\n");
1211 c->SetX(cath, c->GetX(cath)/c->GetCharge(cath));
1213 c->SetX(cath, fSeg2[cath]->GetAnod(fInput->DetElemId(), c->GetX(cath)));
1214 c->SetY(cath, c->GetY(cath)/c->GetCharge(cath));
1216 // apply correction to the coordinate along the anode wire
1222 fSeg2[cath]->GetPadI(fInput->DetElemId(), x, y, fZPlane, ix, iy);
1223 fSeg2[cath]->GetPadC(fInput->DetElemId(), ix, iy, x, y, z);
1224 isec=fSeg2[cath]->Sector(fInput->DetElemId(), ix,iy);
1225 cogCorr = fSeg2[cath]->CorrFunc(fInput->DetElemId(), isec-1);
1230 yOnPad=(c->GetY(cath)-y)/fSeg2[cath]->Dpy(fInput->DetElemId(), isec);
1232 c->SetY(cath, c->GetY(cath)-cogCorr->Eval(yOnPad, 0, 0));
1233 // slat ID from digit
1239 void AliMUONClusterFinderVS::FillCluster(AliMUONRawCluster* c, Int_t cath)
1242 // Completes cluster information starting from list of digits
1252 Float_t xpad, ypad, zpad;
1255 for (Int_t i=0; i<c->GetMultiplicity(cath); i++)
1257 dig = fInput->Digit(cath,c->GetIndex(i,cath));
1259 GetPadC(fInput->DetElemId(),dig->PadX(),dig->PadY(),xpad,ypad, zpad);
1260 AliDebug(1,Form("x %f y %f cx %f cy %f\n",xpad,ypad,c->GetX(0),c->GetY(0)));
1261 dx = xpad - c->GetX(0);
1262 dy = ypad - c->GetY(0);
1263 dr = TMath::Sqrt(dx*dx+dy*dy);
1267 AliDebug(1,Form(" dr %f\n",dr));
1268 Int_t q=dig->Signal();
1269 if (dig->Physics() >= dig->Signal()) {
1271 } else if (dig->Physics() == 0) {
1273 } else c->SetPhysics(i,1);
1274 c->SetPeakSignal(cath,q);
1275 c->SetTrack(0,dig->Hit());
1276 c->SetTrack(1,dig->Track(0));
1277 c->SetTrack(2,dig->Track(1));
1279 AliDebug(1,Form(" c->fTracks[0] %d c->fTracks[1] %d\n",dig->Hit(),
1283 } // loop over digits
1285 // apply correction to the coordinate along the anode wire
1287 c->SetX(cath,fSeg2[cath]->GetAnod(fInput->DetElemId(), c->GetX(cath)));
1290 void AliMUONClusterFinderVS::FindCluster(Int_t i, Int_t j, Int_t cath, AliMUONRawCluster &c){
1294 // Find a super cluster on both cathodes
1297 // Add i,j as element of the cluster
1300 Int_t idx = fHitMap[cath]->GetHitIndex(i,j);
1301 AliMUONDigit* dig = (AliMUONDigit*) fHitMap[cath]->GetHit(i,j);
1302 Int_t q=dig->Signal();
1303 Int_t theX=dig->PadX();
1304 Int_t theY=dig->PadY();
1306 if (q > TMath::Abs(c.GetPeakSignal(0)) && q > TMath::Abs(c.GetPeakSignal(1))) {
1307 c.SetPeakSignal(cath,q);
1308 c.SetTrack(0,dig->Hit());
1309 c.SetTrack(1,dig->Track(0));
1310 c.SetTrack(2,dig->Track(1));
1314 // Make sure that list of digits is ordered
1316 Int_t mu=c.GetMultiplicity(cath);
1317 c.SetIndex(mu, cath, idx);
1319 if (dig->Physics() >= dig->Signal()) {
1321 } else if (dig->Physics() == 0) {
1323 } else c.SetPhysics(mu,1);
1327 for (Int_t ind = mu-1; ind >= 0; ind--) {
1328 Int_t ist=c.GetIndex(ind,cath);
1329 Int_t ql=fInput->Digit(cath, ist)->Signal();
1330 Int_t ix=fInput->Digit(cath, ist)->PadX();
1331 Int_t iy=fInput->Digit(cath, ist)->PadY();
1333 if (q>ql || (q==ql && theX > ix && theY < iy)) {
1334 c.SetIndex(ind, cath, idx);
1335 c.SetIndex(ind+1, cath, ist);
1343 c.SetMultiplicity(cath, c.GetMultiplicity(cath)+1);
1344 if (c.GetMultiplicity(cath) >= 50 ) {
1345 AliDebug(1,Form("FindCluster - multiplicity >50 %d \n",c.GetMultiplicity(0)));
1346 c.SetMultiplicity(cath, 49);
1349 // Prepare center of gravity calculation
1351 fSeg2[cath]->GetPadC(fInput->DetElemId(), i, j, x, y, z);
1354 c.AddCharge(cath,q);
1356 // Flag hit as "taken"
1357 fHitMap[cath]->FlagHit(i,j);
1359 // Now look recursively for all neighbours and pad hit on opposite cathode
1361 // Loop over neighbours
1365 Int_t xList[10], yList[10];
1366 fSeg2[cath]->Neighbours(fInput->DetElemId(), i,j,&nn,xList,yList);
1367 for (Int_t in=0; in<nn; in++) {
1371 if (fHitMap[cath]->TestHit(ix,iy)==kUnused) {
1372 AliDebug(2,Form("\n Neighbours %d %d %d", cath, ix, iy));
1373 FindCluster(ix, iy, cath, c);
1378 Int_t iXopp[50], iYopp[50];
1380 // Neighbours on opposite cathode
1381 // Take into account that several pads can overlap with the present pad
1383 isec=fSeg2[cath]->Sector(fInput->DetElemId(), i,j);
1390 dx = (fSeg2[cath]->Dpx(fInput->DetElemId(), isec))/2.;
1395 dy = (fSeg2[cath]->Dpy(fInput->DetElemId(), isec))/2;
1400 // loop over pad neighbours on opposite cathode
1401 for (fSeg2[iop]->FirstPad(fInput->DetElemId(), x, y, fZPlane, dx, dy);
1402 fSeg2[iop]->MorePads(fInput->DetElemId());
1403 fSeg2[iop]->NextPad(fInput->DetElemId()))
1406 ix = fSeg2[iop]->Ix(); iy = fSeg2[iop]->Iy();
1407 AliDebug(2,Form("\n ix, iy: %f %f %f %d %d %d", x,y,z,ix, iy, fSector));
1408 if (fHitMap[iop]->TestHit(ix,iy)==kUnused){
1411 AliDebug(2,Form("\n Opposite %d %d %d", iop, ix, iy));
1414 } // Loop over pad neighbours
1415 // This had to go outside the loop since recursive calls inside the iterator are not possible
1418 for (jopp=0; jopp<nOpp; jopp++) {
1419 if (fHitMap[iop]->TestHit(iXopp[jopp],iYopp[jopp]) == kUnused)
1420 FindCluster(iXopp[jopp], iYopp[jopp], iop, c);
1425 //_____________________________________________________________________________
1427 void AliMUONClusterFinderVS::FindRawClusters()
1430 // MUON cluster finder from digits -- finds neighbours on both cathodes and
1431 // fills the tree with raw clusters
1435 // Return if no input datad available
1436 if (!fInput->NDigits(0) && !fInput->NDigits(1)) return;
1438 fSeg2[0] = fInput->Segmentation2(0);
1439 fSeg2[1] = fInput->Segmentation2(1);
1441 fHitMap[0] = new AliMUONHitMapA1(fInput->DetElemId(), fSeg2[0], fInput->Digits(0));
1442 fHitMap[1] = new AliMUONHitMapA1(fInput->DetElemId(), fSeg2[1], fInput->Digits(1));
1449 fHitMap[0]->FillHits();
1450 fHitMap[1]->FillHits();
1452 // Outer Loop over Cathodes
1453 for (cath=0; cath<2; cath++) {
1455 for (ndig=0; ndig<fInput->NDigits(cath); ndig++) {
1456 dig = fInput->Digit(cath, ndig);
1457 Int_t padx = dig->PadX();
1458 Int_t pady = dig->PadY();
1459 if (fHitMap[cath]->TestHit(padx,pady)==kUsed ||fHitMap[0]->TestHit(padx,pady)==kEmpty) {
1463 AliDebug(1,Form("\n CATHODE %d CLUSTER %d\n",cath,ncls));
1464 AliMUONRawCluster clus;
1465 clus.SetMultiplicity(0, 0);
1466 clus.SetMultiplicity(1, 0);
1467 clus.SetPeakSignal(cath,dig->Signal());
1468 clus.SetTrack(0, dig->Hit());
1469 clus.SetTrack(1, dig->Track(0));
1470 clus.SetTrack(2, dig->Track(1));
1472 AliDebug(1,Form("idDE %d Padx %d Pady %d", fInput->DetElemId(), padx, pady));
1474 // tag the beginning of cluster list in a raw cluster
1475 clus.SetNcluster(0,-1);
1477 fSeg2[cath]->GetPadC(fInput->DetElemId(), padx, pady, xcu, ycu, fZPlane);
1478 fSector= fSeg2[cath]->Sector(fInput->DetElemId(), padx, pady)/100;
1483 FindCluster(padx,pady,cath,clus);
1484 //^^^^^^^^^^^^^^^^^^^^^^^^
1485 // center of gravity
1486 if (clus.GetX(0)!=0.) clus.SetX(0, clus.GetX(0)/clus.GetCharge(0)); // clus.fX[0] /= clus.fQ[0];
1489 clus.SetX(0,fSeg2[0]->GetAnod(fInput->DetElemId(), clus.GetX(0)));
1490 if (clus.GetY(0)!=0.) clus.SetY(0, clus.GetY(0)/clus.GetCharge(0)); // clus.fY[0] /= clus.fQ[0];
1492 if(clus.GetCharge(1)!=0.) clus.SetX(1, clus.GetX(1)/clus.GetCharge(1)); // clus.fX[1] /= clus.fQ[1];
1495 clus.SetX(1, fSeg2[0]->GetAnod(fInput->DetElemId(),clus.GetX(1)));
1496 if(clus.GetCharge(1)!=0.) clus.SetY(1, clus.GetY(1)/clus.GetCharge(1));// clus.fY[1] /= clus.fQ[1];
1498 clus.SetZ(0, fZPlane);
1499 clus.SetZ(1, fZPlane);
1501 AliDebug(1,Form("\n Cathode 1 multiplicite %d X(CG) %f Y(CG) %f\n",
1502 clus.GetMultiplicity(0),clus.GetX(0),clus.GetY(0)));
1503 AliDebug(1,Form(" Cathode 2 multiplicite %d X(CG) %f Y(CG) %f\n",
1504 clus.GetMultiplicity(1),clus.GetX(1),clus.GetY(1)));
1505 // Analyse cluster and decluster if necessary
1508 clus.SetNcluster(1,fNRawClusters);
1509 clus.SetClusterType(clus.PhysicsContribution());
1516 // reset Cluster object
1517 { // begin local scope
1518 for (int k=0;k<clus.GetMultiplicity(0);k++) clus.SetIndex(k, 0, 0);
1519 } // end local scope
1521 { // begin local scope
1522 for (int k=0;k<clus.GetMultiplicity(1);k++) clus.SetIndex(k, 1, 0);
1523 } // end local scope
1525 clus.SetMultiplicity(0,0);
1526 clus.SetMultiplicity(1,0);
1530 } // end loop cathodes
1535 Float_t AliMUONClusterFinderVS::SingleMathiesonFit(AliMUONRawCluster *c, Int_t cath)
1537 // Performs a single Mathieson fit on one cathode
1539 Double_t arglist[20];
1541 AliMUONClusterInput& clusterInput = *(AliMUONClusterInput::Instance());
1543 clusterInput.Fitter()->SetFCN(fcnS1);
1544 clusterInput.Fitter()->mninit(2,10,7);
1545 clusterInput.Fitter()->SetPrintLevel(-1 + AliLog::GetGlobalDebugLevel());
1547 clusterInput.Fitter()->mnexcm("SET NOW", arglist, 0, ierflag);
1548 // Set starting values
1549 static Double_t vstart[2];
1550 vstart[0]=c->GetX(1);
1551 vstart[1]=c->GetY(0);
1554 // lower and upper limits
1555 static Double_t lower[2], upper[2];
1557 fSeg2[cath]->GetPadI(fInput->DetElemId(), c->GetX(cath), c->GetY(cath), fZPlane, ix, iy);
1558 isec=fSeg2[cath]->Sector(fInput->DetElemId(), ix, iy);
1560 lower[0]=vstart[0]-fSeg2[cath]->Dpx(fInput->DetElemId(), isec)/2;
1561 lower[1]=vstart[1]-fSeg2[cath]->Dpy(fInput->DetElemId(), isec)/2;
1563 upper[0]=lower[0]+fSeg2[cath]->Dpx(fInput->DetElemId(), isec);
1564 upper[1]=lower[1]+fSeg2[cath]->Dpy(fInput->DetElemId(), isec);
1568 static Double_t step[2]={0.0005, 0.0005};
1570 clusterInput.Fitter()->mnparm(0,"x1",vstart[0],step[0],lower[0],upper[0],ierflag);
1571 clusterInput.Fitter()->mnparm(1,"y1",vstart[1],step[1],lower[1],upper[1],ierflag);
1572 // ready for minimisation
1576 clusterInput.Fitter()->mnexcm("SET NOGR", arglist, 0, ierflag);
1577 clusterInput.Fitter()->mnexcm("MIGRAD", arglist, 0, ierflag);
1578 // clusterInput.Fitter()->mnexcm("EXIT" , arglist, 0, ierflag);
1579 Double_t fmin, fedm, errdef;
1580 Int_t npari, nparx, istat;
1582 clusterInput.Fitter()->mnstat(fmin, fedm, errdef, npari, nparx, istat);
1586 // Get fitted parameters
1587 Double_t xrec, yrec;
1589 Double_t epxz, b1, b2;
1591 clusterInput.Fitter()->mnpout(0, chname, xrec, epxz, b1, b2, ierflg);
1592 clusterInput.Fitter()->mnpout(1, chname, yrec, epxz, b1, b2, ierflg);
1598 Float_t AliMUONClusterFinderVS::CombiSingleMathiesonFit(AliMUONRawCluster * /*c*/)
1600 // Perform combined Mathieson fit on both cathode planes
1602 Double_t arglist[20];
1604 AliMUONClusterInput& clusterInput = *(AliMUONClusterInput::Instance());
1605 clusterInput.Fitter()->SetFCN(fcnCombiS1);
1606 clusterInput.Fitter()->mninit(2,10,7);
1607 clusterInput.Fitter()->SetPrintLevel(-1 + AliLog::GetGlobalDebugLevel());
1609 clusterInput.Fitter()->mnexcm("SET NOW", arglist, 0, ierflag);
1610 static Double_t vstart[2];
1611 vstart[0]=fXInit[0];
1612 vstart[1]=fYInit[0];
1615 // lower and upper limits
1616 static Float_t lower[2], upper[2];
1620 fSeg2[0]->GetPadI(fInput->DetElemId(), fXInit[0], fYInit[0], fZPlane, ix, iy);
1621 isec=fSeg2[0]->Sector(fInput->DetElemId(), ix, iy);
1622 dpy=fSeg2[0]->Dpy(fInput->DetElemId(), isec);
1623 fSeg2[1]->GetPadI(fInput->DetElemId(), fXInit[0], fYInit[0], fZPlane, ix, iy);
1624 isec=fSeg2[1]->Sector(fInput->DetElemId(), ix, iy);
1625 dpx=fSeg2[1]->Dpx(fInput->DetElemId(), isec);
1628 Float_t xdum, ydum, zdum;
1630 // Find save upper and lower limits
1633 for (fSeg2[1]->FirstPad(fInput->DetElemId(),fXInit[0], fYInit[0], fZPlane, dpx, 0.);
1634 fSeg2[1]->MorePads(fInput->DetElemId());
1635 fSeg2[1]->NextPad(fInput->DetElemId()))
1637 ix=fSeg2[1]->Ix(); iy=fSeg2[1]->Iy();
1638 fSeg2[1]->GetPadC(fInput->DetElemId(), ix,iy, upper[0], ydum, zdum);
1639 if (icount ==0) lower[0]=upper[0];
1643 if (lower[0]>upper[0]) {xdum=lower[0]; lower[0]=upper[0]; upper[0]=xdum;}
1646 AliDebug(1,Form("\n single y %f %f", fXInit[0], fYInit[0]));
1648 for (fSeg2[0]->FirstPad(fInput->DetElemId(), fXInit[0], fYInit[0], fZPlane, 0., dpy);
1649 fSeg2[0]->MorePads(fInput->DetElemId());
1650 fSeg2[0]->NextPad(fInput->DetElemId()))
1652 ix=fSeg2[0]->Ix(); iy=fSeg2[0]->Iy();
1653 fSeg2[0]->GetPadC(fInput->DetElemId(), ix,iy,xdum,upper[1],zdum);
1654 if (icount ==0) lower[1]=upper[1];
1656 AliDebug(1,Form("\n upper lower %d %f %f", icount, upper[1], lower[1]));
1659 if (lower[1]>upper[1]) {xdum=lower[1]; lower[1]=upper[1]; upper[1]=xdum;}
1662 static Double_t step[2]={0.00001, 0.0001};
1664 clusterInput.Fitter()->mnparm(0,"x1",vstart[0],step[0],lower[0],upper[0],ierflag);
1665 clusterInput.Fitter()->mnparm(1,"y1",vstart[1],step[1],lower[1],upper[1],ierflag);
1666 // ready for minimisation
1670 clusterInput.Fitter()->mnexcm("SET NOGR", arglist, 0, ierflag);
1671 clusterInput.Fitter()->mnexcm("MIGRAD", arglist, 0, ierflag);
1672 // clusterInput.Fitter()->mnexcm("EXIT" , arglist, 0, ierflag);
1673 Double_t fmin, fedm, errdef;
1674 Int_t npari, nparx, istat;
1676 clusterInput.Fitter()->mnstat(fmin, fedm, errdef, npari, nparx, istat);
1680 // Get fitted parameters
1681 Double_t xrec, yrec;
1683 Double_t epxz, b1, b2;
1685 clusterInput.Fitter()->mnpout(0, chname, xrec, epxz, b1, b2, ierflg);
1686 clusterInput.Fitter()->mnpout(1, chname, yrec, epxz, b1, b2, ierflg);
1692 Bool_t AliMUONClusterFinderVS::DoubleMathiesonFit(AliMUONRawCluster * /*c*/, Int_t cath)
1694 // Performs a double Mathieson fit on one cathode
1698 // Initialise global variables for fit
1699 Double_t arglist[20];
1701 AliMUONClusterInput& clusterInput = *(AliMUONClusterInput::Instance());
1702 clusterInput.Fitter()->SetFCN(fcnS2);
1703 clusterInput.Fitter()->mninit(5,10,7);
1704 clusterInput.Fitter()->SetPrintLevel(-1 + AliLog::GetGlobalDebugLevel());
1706 clusterInput.Fitter()->mnexcm("SET NOW", arglist, 0, ierflag);
1707 // Set starting values
1708 static Double_t vstart[5];
1709 vstart[0]=fX[fIndLocal[0][cath]][cath];
1710 vstart[1]=fY[fIndLocal[0][cath]][cath];
1711 vstart[2]=fX[fIndLocal[1][cath]][cath];
1712 vstart[3]=fY[fIndLocal[1][cath]][cath];
1713 vstart[4]=Float_t(fQ[fIndLocal[0][cath]][cath])/
1714 Float_t(fQ[fIndLocal[0][cath]][cath]+fQ[fIndLocal[1][cath]][cath]);
1715 // lower and upper limits
1716 static Float_t lower[5], upper[5];
1719 isec=fSeg2[cath]->Sector(fInput->DetElemId(),fIx[fIndLocal[0][cath]][cath],
1720 fIy[fIndLocal[0][cath]][cath]);
1721 lower[0]=vstart[0]-fSeg2[cath]->Dpx(fInput->DetElemId(),isec);
1722 lower[1]=vstart[1]-fSeg2[cath]->Dpy(fInput->DetElemId(),isec);
1724 upper[0]=lower[0]+2.*fSeg2[cath]->Dpx(fInput->DetElemId(),isec);
1725 upper[1]=lower[1]+2.*fSeg2[cath]->Dpy(fInput->DetElemId(),isec);
1727 isec=fSeg2[cath]->Sector(fInput->DetElemId(),fIx[fIndLocal[1][cath]][cath],
1728 fIy[fIndLocal[1][cath]][cath]);
1729 lower[2]=vstart[2]-fSeg2[cath]->Dpx(fInput->DetElemId(),isec)/2;
1730 lower[3]=vstart[3]-fSeg2[cath]->Dpy(fInput->DetElemId(),isec)/2;
1732 upper[2]=lower[2]+fSeg2[cath]->Dpx(fInput->DetElemId(),isec);
1733 upper[1]=lower[1]+2.*fSeg2[cath]->Dpy(fInput->DetElemId(),isec);
1740 static Double_t step[5]={0.0005, 0.0005, 0.0005, 0.0005, 0.0001};
1742 clusterInput.Fitter()->mnparm(0,"x1",vstart[0],step[0],lower[0],upper[0],ierflag);
1743 clusterInput.Fitter()->mnparm(1,"y1",vstart[1],step[1],lower[1],upper[1],ierflag);
1744 clusterInput.Fitter()->mnparm(2,"x2",vstart[2],step[2],lower[2],upper[2],ierflag);
1745 clusterInput.Fitter()->mnparm(3,"y2",vstart[3],step[3],lower[3],upper[3],ierflag);
1746 clusterInput.Fitter()->mnparm(4,"a0",vstart[4],step[4],lower[4],upper[4],ierflag);
1747 // ready for minimisation
1751 clusterInput.Fitter()->mnexcm("SET NOGR", arglist, 0, ierflag);
1752 clusterInput.Fitter()->mnexcm("MIGRAD", arglist, 0, ierflag);
1753 // clusterInput.Fitter()->mnexcm("EXIT" , arglist, 0, ierflag);
1754 // Get fitted parameters
1755 Double_t xrec[2], yrec[2], qfrac;
1757 Double_t epxz, b1, b2;
1759 clusterInput.Fitter()->mnpout(0, chname, xrec[0], epxz, b1, b2, ierflg);
1760 clusterInput.Fitter()->mnpout(1, chname, yrec[0], epxz, b1, b2, ierflg);
1761 clusterInput.Fitter()->mnpout(2, chname, xrec[1], epxz, b1, b2, ierflg);
1762 clusterInput.Fitter()->mnpout(3, chname, yrec[1], epxz, b1, b2, ierflg);
1763 clusterInput.Fitter()->mnpout(4, chname, qfrac, epxz, b1, b2, ierflg);
1765 Double_t fmin, fedm, errdef;
1766 Int_t npari, nparx, istat;
1768 clusterInput.Fitter()->mnstat(fmin, fedm, errdef, npari, nparx, istat);
1773 Float_t AliMUONClusterFinderVS::CombiDoubleMathiesonFit(AliMUONRawCluster * /*c*/)
1776 // Perform combined double Mathieson fit on both cathode planes
1778 Double_t arglist[20];
1780 AliMUONClusterInput& clusterInput = *(AliMUONClusterInput::Instance());
1781 clusterInput.Fitter()->SetFCN(fcnCombiS2);
1782 clusterInput.Fitter()->mninit(6,10,7);
1783 clusterInput.Fitter()->SetPrintLevel(-1 + AliLog::GetGlobalDebugLevel());
1785 clusterInput.Fitter()->mnexcm("SET NOW", arglist, 0, ierflag);
1786 // Set starting values
1787 static Double_t vstart[6];
1788 vstart[0]=fXInit[0];
1789 vstart[1]=fYInit[0];
1790 vstart[2]=fXInit[1];
1791 vstart[3]=fYInit[1];
1792 vstart[4]=fQrInit[0];
1793 vstart[5]=fQrInit[1];
1794 // lower and upper limits
1795 static Float_t lower[6], upper[6];
1799 fSeg2[1]->GetPadI(fInput->DetElemId(),fXInit[0], fYInit[0], fZPlane, ix, iy);
1800 isec=fSeg2[1]->Sector(fInput->DetElemId(),ix, iy);
1801 dpx=fSeg2[1]->Dpx(fInput->DetElemId(), isec);
1803 fSeg2[0]->GetPadI(fInput->DetElemId(), fXInit[0], fYInit[0], fZPlane, ix, iy);
1804 isec=fSeg2[0]->Sector(fInput->DetElemId(), ix, iy);
1805 dpy=fSeg2[0]->Dpy(fInput->DetElemId(), isec);
1810 Float_t xdum, ydum, zdum;
1811 AliDebug(1,Form("\n Cluster Finder: %f %f %f %f ", fXInit[0], fXInit[1],fYInit[0], fYInit[1] ));
1813 // Find save upper and lower limits
1816 for (fSeg2[1]->FirstPad(fInput->DetElemId(),fXInit[0], fYInit[0], fZPlane, dpx, 0.);
1817 fSeg2[1]->MorePads(fInput->DetElemId());
1818 fSeg2[1]->NextPad(fInput->DetElemId()))
1820 ix=fSeg2[1]->Ix(); iy=fSeg2[1]->Iy();
1821 // if (fHitMap[1]->TestHit(ix, iy) == kEmpty) continue;
1822 fSeg2[1]->GetPadC(fInput->DetElemId(),ix,iy,upper[0],ydum,zdum);
1823 if (icount ==0) lower[0]=upper[0];
1826 if (lower[0]>upper[0]) {xdum=lower[0]; lower[0]=upper[0]; upper[0]=xdum;}
1827 // vstart[0] = 0.5*(lower[0]+upper[0]);
1832 for (fSeg2[0]->FirstPad(fInput->DetElemId(),fXInit[0], fYInit[0], fZPlane, 0., dpy);
1833 fSeg2[0]->MorePads(fInput->DetElemId());
1834 fSeg2[0]->NextPad(fInput->DetElemId()))
1836 ix=fSeg2[0]->Ix(); iy=fSeg2[0]->Iy();
1837 // if (fHitMap[0]->TestHit(ix, iy) == kEmpty) continue;
1838 fSeg2[0]->GetPadC(fInput->DetElemId(),ix,iy,xdum,upper[1],zdum);
1839 if (icount ==0) lower[1]=upper[1];
1843 if (lower[1]>upper[1]) {xdum=lower[1]; lower[1]=upper[1]; upper[1]=xdum;}
1844 // vstart[1] = 0.5*(lower[1]+upper[1]);
1847 fSeg2[1]->GetPadI(fInput->DetElemId(),fXInit[1], fYInit[1], fZPlane, ix, iy);
1848 isec=fSeg2[1]->Sector(fInput->DetElemId(),ix, iy);
1849 dpx=fSeg2[1]->Dpx(fInput->DetElemId(),isec);
1850 fSeg2[0]->GetPadI(fInput->DetElemId(),fXInit[1], fYInit[1], fZPlane, ix, iy);
1851 isec=fSeg2[0]->Sector(fInput->DetElemId(),ix, iy);
1852 dpy=fSeg2[0]->Dpy(fInput->DetElemId(),isec);
1855 // Find save upper and lower limits
1859 for (fSeg2[1]->FirstPad(fInput->DetElemId(),fXInit[1], fYInit[1], fZPlane, dpx, 0);
1860 fSeg2[1]->MorePads(fInput->DetElemId());
1861 fSeg2[1]->NextPad(fInput->DetElemId()))
1863 ix=fSeg2[1]->Ix(); iy=fSeg2[1]->Iy();
1864 // if (fHitMap[1]->TestHit(ix, iy) == kEmpty) continue;
1865 fSeg2[1]->GetPadC(fInput->DetElemId(),ix,iy,upper[2],ydum,zdum);
1866 if (icount ==0) lower[2]=upper[2];
1869 if (lower[2]>upper[2]) {xdum=lower[2]; lower[2]=upper[2]; upper[2]=xdum;}
1870 // vstart[2] = 0.5*(lower[2]+upper[2]);
1874 for (fSeg2[0]->FirstPad(fInput->DetElemId(),fXInit[1], fYInit[1], fZPlane, 0, dpy);
1875 fSeg2[0]-> MorePads(fInput->DetElemId());
1876 fSeg2[0]->NextPad(fInput->DetElemId()))
1878 ix=fSeg2[0]->Ix(); iy=fSeg2[0]->Iy();
1879 // if (fHitMap[0]->TestHit(ix, iy) != kEmpty) continue;
1881 fSeg2[0]->GetPadC(fInput->DetElemId(),ix,iy,xdum,upper[3],zdum);
1882 if (icount ==0) lower[3]=upper[3];
1886 if (lower[3]>upper[3]) {xdum=lower[3]; lower[3]=upper[3]; upper[3]=xdum;}
1894 static Double_t step[6]={0.0005, 0.0005, 0.0005, 0.0005, 0.001, 0.001};
1895 clusterInput.Fitter()->mnparm(0,"x1",vstart[0],step[0],lower[0],upper[0],ierflag);
1896 clusterInput.Fitter()->mnparm(1,"y1",vstart[1],step[1],lower[1],upper[1],ierflag);
1897 clusterInput.Fitter()->mnparm(2,"x2",vstart[2],step[2],lower[2],upper[2],ierflag);
1898 clusterInput.Fitter()->mnparm(3,"y2",vstart[3],step[3],lower[3],upper[3],ierflag);
1899 clusterInput.Fitter()->mnparm(4,"a0",vstart[4],step[4],lower[4],upper[4],ierflag);
1900 clusterInput.Fitter()->mnparm(5,"a1",vstart[5],step[5],lower[5],upper[5],ierflag);
1901 // ready for minimisation
1905 clusterInput.Fitter()->mnexcm("SET NOGR", arglist, 0, ierflag);
1906 clusterInput.Fitter()->mnexcm("MIGRAD", arglist, 0, ierflag);
1907 // clusterInput.Fitter()->mnexcm("EXIT" , arglist, 0, ierflag);
1908 // Get fitted parameters
1910 Double_t epxz, b1, b2;
1912 clusterInput.Fitter()->mnpout(0, chname, fXFit[0], epxz, b1, b2, ierflg);
1913 clusterInput.Fitter()->mnpout(1, chname, fYFit[0], epxz, b1, b2, ierflg);
1914 clusterInput.Fitter()->mnpout(2, chname, fXFit[1], epxz, b1, b2, ierflg);
1915 clusterInput.Fitter()->mnpout(3, chname, fYFit[1], epxz, b1, b2, ierflg);
1916 clusterInput.Fitter()->mnpout(4, chname, fQrFit[0], epxz, b1, b2, ierflg);
1917 clusterInput.Fitter()->mnpout(5, chname, fQrFit[1], epxz, b1, b2, ierflg);
1919 Double_t fmin, fedm, errdef;
1920 Int_t npari, nparx, istat;
1922 clusterInput.Fitter()->mnstat(fmin, fedm, errdef, npari, nparx, istat);
1930 void AliMUONClusterFinderVS::Split(AliMUONRawCluster* c)
1933 // One cluster for each maximum
1936 AliMUONClusterInput& clusterInput = *(AliMUONClusterInput::Instance());
1937 for (j=0; j<2; j++) {
1938 AliMUONRawCluster cnew;
1939 cnew.SetGhost(c->GetGhost());
1940 for (cath=0; cath<2; cath++) {
1941 cnew.SetChi2(cath,fChi2[0]);
1942 // ?? why not cnew.fChi2[cath]=fChi2[cath];
1945 cnew.SetNcluster(0,-1);
1946 cnew.SetNcluster(1,fNRawClusters);
1948 cnew.SetNcluster(0,fNPeaks);
1949 cnew.SetNcluster(1,0);
1951 cnew.SetMultiplicity(cath,0);
1952 cnew.SetX(cath, Float_t(fXFit[j]));
1953 cnew.SetY(cath, Float_t(fYFit[j]));
1954 cnew.SetZ(cath, fZPlane);
1956 cnew.SetCharge(cath, Int_t(clusterInput.TotalCharge(cath)*fQrFit[cath]));
1958 cnew.SetCharge(cath, Int_t(clusterInput.TotalCharge(cath)*(1-fQrFit[cath])));
1960 fSeg2[cath]->SetHit(fInput->DetElemId(), fXFit[j],fYFit[j],fZPlane);
1962 for (i=0; i<fMul[cath]; i++) {
1964 cnew.SetIndex(cnew.GetMultiplicity(cath), cath, c->GetIndex(i,cath));
1966 fSeg2[cath]->SetPad(fInput->DetElemId(),fIx[i][cath], fIy[i][cath]);
1967 q1 = fInput->Mathieson()->IntXY(fInput->DetElemId(),fSeg2[cath]);
1969 cnew.SetContrib(i, cath, q1*Float_t(cnew.GetCharge(cath))/Float_t(fQ[i][cath]));
1970 cnew.SetMultiplicity(cath, cnew.GetMultiplicity(cath)+1 );
1972 FillCluster(&cnew,0,cath);
1974 cnew.SetClusterType(cnew.PhysicsContribution());
1975 if (cnew.GetCharge(0)>0 && cnew.GetCharge(1)>0) AddRawCluster(cnew);
1979 void AliMUONClusterFinderVS::AddRawCluster(AliMUONRawCluster& c)
1982 // Add a raw cluster copy to the list
1984 // AliMUON *pMUON=(AliMUON*)gAlice->GetModule("MUON");
1985 // pMUON->GetMUONData()->AddRawCluster(fInput->Chamber(),c);
1988 // Setting detection element in raw cluster for alignment
1990 c.SetDetElemId(fInput->DetElemId());
1992 TClonesArray &lrawcl = *fRawClusters;
1993 new(lrawcl[fNRawClusters++]) AliMUONRawCluster(c);
1994 AliDebug(1,Form("\nfNRawClusters %d\n",fNRawClusters));
1997 AliMUONClusterFinderVS& AliMUONClusterFinderVS
1998 ::operator = (const AliMUONClusterFinderVS& rhs)
2000 // Protected assignement operator
2002 if (this == &rhs) return *this;
2004 AliFatal("Not implemented.");
2010 // Minimisation functions
2012 void fcnS1(Int_t & /*npar*/, Double_t * /*gin*/, Double_t &f, Double_t *par, Int_t /*iflag*/)
2014 AliMUONClusterInput& clusterInput = *(AliMUONClusterInput::Instance());
2021 for (i=0; i<clusterInput.Nmul(0); i++) {
2022 Float_t q0=clusterInput.Charge(i,0);
2023 Float_t q1=clusterInput.DiscrChargeS1(i,par);
2032 void fcnCombiS1(Int_t & /*npar*/, Double_t * /*gin*/, Double_t &f, Double_t *par, Int_t /*iflag*/)
2034 AliMUONClusterInput& clusterInput = *(AliMUONClusterInput::Instance());
2041 for (cath=0; cath<2; cath++) {
2042 for (i=0; i<clusterInput.Nmul(cath); i++) {
2043 Float_t q0=clusterInput.Charge(i,cath);
2044 Float_t q1=clusterInput.DiscrChargeCombiS1(i,par,cath);
2055 void fcnS2(Int_t & /*npar*/, Double_t * /*gin*/, Double_t &f, Double_t *par, Int_t /*iflag*/)
2057 AliMUONClusterInput& clusterInput = *(AliMUONClusterInput::Instance());
2064 for (i=0; i<clusterInput.Nmul(0); i++) {
2066 Float_t q0=clusterInput.Charge(i,0);
2067 Float_t q1=clusterInput.DiscrChargeS2(i,par);
2077 void fcnCombiS2(Int_t & /*npar*/, Double_t * /*gin*/, Double_t &f, Double_t *par, Int_t /*iflag*/)
2079 AliMUONClusterInput& clusterInput = *(AliMUONClusterInput::Instance());
2085 for (cath=0; cath<2; cath++) {
2086 for (i=0; i<clusterInput.Nmul(cath); i++) {
2087 Float_t q0=clusterInput.Charge(i,cath);
2088 Float_t q1=clusterInput.DiscrChargeCombiS2(i,par,cath);