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 "AliSegmentation.h"
22 #include "AliMUONResponse.h"
23 #include "AliMUONClusterInput.h"
24 #include "AliMUONHitMapA1.h"
33 #include <TPostScript.h>
38 #include <Riostream.h>
40 //_____________________________________________________________________
41 // This function is minimized in the double-Mathieson fit
42 void fcnS2(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag);
43 void fcnS1(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag);
44 void fcnCombiS1(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag);
45 void fcnCombiS2(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag);
47 ClassImp(AliMUONClusterFinderVS)
49 AliMUONClusterFinderVS::AliMUONClusterFinderVS()
51 // Default constructor
52 fInput=AliMUONClusterInput::Instance();
55 fTrack[0]=fTrack[1]=-1;
56 fDebugLevel = 0; // make silent default
57 fGhostChi2Cut = 1e6; // nothing done by default
60 for(Int_t i=0; i<100; i++) {
61 for (Int_t j=0; j<2; j++) {
65 fRawClusters = new TClonesArray("AliMUONRawCluster",1000);
70 //____________________________________________________________________________
71 AliMUONClusterFinderVS::~AliMUONClusterFinderVS()
73 // Reset tracks information
76 fRawClusters->Delete();
81 AliMUONClusterFinderVS::AliMUONClusterFinderVS(const AliMUONClusterFinderVS & clusterFinder):TObject(clusterFinder)
83 // Dummy copy Constructor
86 //____________________________________________________________________________
87 void AliMUONClusterFinderVS::ResetRawClusters()
89 // Reset tracks information
91 if (fRawClusters) fRawClusters->Clear();
93 //____________________________________________________________________________
94 void AliMUONClusterFinderVS::Decluster(AliMUONRawCluster *cluster)
96 // Decluster by local maxima
97 SplitByLocalMaxima(cluster);
99 //____________________________________________________________________________
100 void AliMUONClusterFinderVS::SplitByLocalMaxima(AliMUONRawCluster *c)
102 // Split complex cluster by local maxima
105 fInput->SetCluster(c);
107 fMul[0]=c->fMultiplicity[0];
108 fMul[1]=c->fMultiplicity[1];
111 // dump digit information into arrays
116 for (cath=0; cath<2; cath++) {
118 for (i=0; i<fMul[cath]; i++)
121 fDig[i][cath]=fInput->Digit(cath, c->fIndexMap[i][cath]);
123 fIx[i][cath]= fDig[i][cath]->PadX();
124 fIy[i][cath]= fDig[i][cath]->PadY();
126 fQ[i][cath] = fDig[i][cath]->Signal();
127 // pad centre coordinates
129 GetPadC(fIx[i][cath], fIy[i][cath], fX[i][cath], fY[i][cath], fZ[i][cath]);
130 } // loop over cluster digits
131 } // loop over cathodes
137 // Initialise and perform mathieson fits
138 Float_t chi2, oldchi2;
139 // ++++++++++++++++++*************+++++++++++++++++++++
140 // (1) No more than one local maximum per cathode plane
141 // +++++++++++++++++++++++++++++++*************++++++++
142 if ((fNLocal[0]==1 && (fNLocal[1]==0 || fNLocal[1]==1)) ||
143 (fNLocal[0]==0 && fNLocal[1]==1)) {
144 // Perform combined single Mathieson fit
145 // Initial values for coordinates (x,y)
147 // One local maximum on cathodes 1 and 2 (X->cathode 2, Y->cathode 1)
148 if (fNLocal[0]==1 && fNLocal[1]==1) {
151 // One local maximum on cathode 1 (X,Y->cathode 1)
152 } else if (fNLocal[0]==1) {
155 // One local maximum on cathode 2 (X,Y->cathode 2)
161 fprintf(stderr,"\n cas (1) CombiSingleMathiesonFit(c)\n");
162 chi2=CombiSingleMathiesonFit(c);
163 // Int_t ndf = fgNbins[0]+fgNbins[1]-2;
164 // Float_t prob = TMath::Prob(Double_t(chi2),ndf);
165 // prob1->Fill(prob);
166 // chi2_1->Fill(chi2);
169 fprintf(stderr," chi2 %f ",chi2);
179 c->fX[0]=fSeg[0]->GetAnod(c->fX[0]);
180 c->fX[1]=fSeg[1]->GetAnod(c->fX[1]);
182 // If reasonable chi^2 add result to the list of rawclusters
185 // If not try combined double Mathieson Fit
188 fprintf(stderr," MAUVAIS CHI2 !!!\n");
189 if (fNLocal[0]==1 && fNLocal[1]==1) {
190 fXInit[0]=fX[fIndLocal[0][1]][1];
191 fYInit[0]=fY[fIndLocal[0][0]][0];
192 fXInit[1]=fX[fIndLocal[0][1]][1];
193 fYInit[1]=fY[fIndLocal[0][0]][0];
194 } else if (fNLocal[0]==1) {
195 fXInit[0]=fX[fIndLocal[0][0]][0];
196 fYInit[0]=fY[fIndLocal[0][0]][0];
197 fXInit[1]=fX[fIndLocal[0][0]][0];
198 fYInit[1]=fY[fIndLocal[0][0]][0];
200 fXInit[0]=fX[fIndLocal[0][1]][1];
201 fYInit[0]=fY[fIndLocal[0][1]][1];
202 fXInit[1]=fX[fIndLocal[0][1]][1];
203 fYInit[1]=fY[fIndLocal[0][1]][1];
206 // Initial value for charge ratios
210 fprintf(stderr,"\n cas (1) CombiDoubleMathiesonFit(c)\n");
211 chi2=CombiDoubleMathiesonFit(c);
212 // Int_t ndf = fgNbins[0]+fgNbins[1]-6;
213 // Float_t prob = TMath::Prob(chi2,ndf);
214 // prob2->Fill(prob);
215 // chi2_2->Fill(chi2);
217 // Was this any better ??
219 fprintf(stderr," Old and new chi2 %f %f ", oldchi2, chi2);
220 if (fFitStat!=0 && chi2>0 && (2.*chi2 < oldchi2)) {
222 fprintf(stderr," Split\n");
223 // Split cluster into two according to fit result
227 fprintf(stderr," Don't Split\n");
233 // +++++++++++++++++++++++++++++++++++++++
234 // (2) Two local maxima per cathode plane
235 // +++++++++++++++++++++++++++++++++++++++
236 } else if (fNLocal[0]==2 && fNLocal[1]==2) {
238 // Let's look for ghosts first
240 Float_t xm[4][2], ym[4][2];
241 Float_t dpx, dpy, dx, dy;
242 Int_t ixm[4][2], iym[4][2];
243 Int_t isec, im1, im2, ico;
245 // Form the 2x2 combinations
246 // 0-0, 0-1, 1-0, 1-1
248 for (im1=0; im1<2; im1++) {
249 for (im2=0; im2<2; im2++) {
250 xm[ico][0]=fX[fIndLocal[im1][0]][0];
251 ym[ico][0]=fY[fIndLocal[im1][0]][0];
252 xm[ico][1]=fX[fIndLocal[im2][1]][1];
253 ym[ico][1]=fY[fIndLocal[im2][1]][1];
255 ixm[ico][0]=fIx[fIndLocal[im1][0]][0];
256 iym[ico][0]=fIy[fIndLocal[im1][0]][0];
257 ixm[ico][1]=fIx[fIndLocal[im2][1]][1];
258 iym[ico][1]=fIy[fIndLocal[im2][1]][1];
262 // ico = 0 : first local maximum on cathodes 1 and 2
263 // ico = 1 : fisrt local maximum on cathode 1 and second on cathode 2
264 // ico = 2 : second local maximum on cathode 1 and first on cathode 1
265 // ico = 3 : second local maximum on cathodes 1 and 2
267 // Analyse the combinations and keep those that are possible !
268 // For each combination check consistency in x and y
271 Float_t dr[4] = {1.e4, 1.e4, 1.e4, 1.e4};
274 // In case of staggering maxima are displaced by exactly half the pad-size in y.
275 // We have to take into account the numerical precision in the consistency check;
278 for (ico=0; ico<4; ico++) {
279 accepted[ico]=kFALSE;
280 // cathode one: x-coordinate
281 isec=fSeg[0]->Sector(ixm[ico][0], iym[ico][0]);
282 dpx=fSeg[0]->Dpx(isec)/2.;
283 dx=TMath::Abs(xm[ico][0]-xm[ico][1]);
284 // cathode two: y-coordinate
285 isec=fSeg[1]->Sector(ixm[ico][1], iym[ico][1]);
286 dpy=fSeg[1]->Dpy(isec)/2.;
287 dy=TMath::Abs(ym[ico][0]-ym[ico][1]);
289 printf("\n %i %f %f %f %f %f %f \n", ico, ym[ico][0], ym[ico][1], dy, dpy, dx, dpx );
290 if ((dx <= dpx) && (dy <= dpy+eps)) {
293 dr[ico] = TMath::Sqrt(dx*dx+dy*dy);
297 accepted[ico]=kFALSE;
301 printf("\n iacc= %d:\n", iacc);
303 if (accepted[0] && accepted[1]) {
304 if (dr[0] >= dr[1]) {
311 if (accepted[2] && accepted[3]) {
312 if (dr[2] >= dr[3]) {
319 // eliminate one candidate
323 for (ico=0; ico<4; ico++) {
324 if (accepted[ico] && dr[ico] > drmax) {
330 accepted[icobad] = kFALSE;
337 printf("\n iacc= %d:\n", iacc);
339 fprintf(stderr,"\n iacc=2: No problem ! \n");
340 } else if (iacc==4) {
341 fprintf(stderr,"\n iacc=4: Ok, but ghost problem !!! \n");
342 } else if (iacc==0) {
343 fprintf(stderr,"\n iacc=0: I don't know what to do with this !!!!!!!!! \n");
347 // Initial value for charge ratios
348 fQrInit[0]=Float_t(fQ[fIndLocal[0][0]][0])/
349 Float_t(fQ[fIndLocal[0][0]][0]+fQ[fIndLocal[1][0]][0]);
350 fQrInit[1]=Float_t(fQ[fIndLocal[0][1]][1])/
351 Float_t(fQ[fIndLocal[0][1]][1]+fQ[fIndLocal[1][1]][1]);
353 // ******* iacc = 0 *******
354 // No combinations found between the 2 cathodes
355 // We keep the center of gravity of the cluster
360 // ******* iacc = 1 *******
361 // Only one combination found between the 2 cathodes
363 // Initial values for the 2 maxima (x,y)
365 // 1 maximum is initialised with the maximum of the combination found (X->cathode 2, Y->cathode 1)
366 // 1 maximum is initialised with the other maximum of the first cathode
368 fprintf(stderr,"ico=0\n");
373 } else if (accepted[1]){
374 fprintf(stderr,"ico=1\n");
379 } else if (accepted[2]){
380 fprintf(stderr,"ico=2\n");
385 } else if (accepted[3]){
386 fprintf(stderr,"ico=3\n");
393 fprintf(stderr,"\n cas (2) CombiDoubleMathiesonFit(c)\n");
394 chi2=CombiDoubleMathiesonFit(c);
395 // Int_t ndf = fgNbins[0]+fgNbins[1]-6;
396 // Float_t prob = TMath::Prob(chi2,ndf);
397 // prob2->Fill(prob);
398 // chi2_2->Fill(chi2);
400 fprintf(stderr," chi2 %f\n",chi2);
402 // If reasonable chi^2 add result to the list of rawclusters
407 // 1 maximum is initialised with the maximum of the combination found (X->cathode 2, Y->cathode 1)
408 // 1 maximum is initialised with the other maximum of the second cathode
410 fprintf(stderr,"ico=0\n");
415 } else if (accepted[1]){
416 fprintf(stderr,"ico=1\n");
421 } else if (accepted[2]){
422 fprintf(stderr,"ico=2\n");
427 } else if (accepted[3]){
428 fprintf(stderr,"ico=3\n");
435 fprintf(stderr,"\n cas (2) CombiDoubleMathiesonFit(c)\n");
436 chi2=CombiDoubleMathiesonFit(c);
437 // Int_t ndf = fgNbins[0]+fgNbins[1]-6;
438 // Float_t prob = TMath::Prob(chi2,ndf);
439 // prob2->Fill(prob);
440 // chi2_2->Fill(chi2);
442 fprintf(stderr," chi2 %f\n",chi2);
444 // If reasonable chi^2 add result to the list of rawclusters
448 //We keep only the combination found (X->cathode 2, Y->cathode 1)
449 for (Int_t ico=0; ico<2; ico++) {
451 AliMUONRawCluster cnew;
453 for (cath=0; cath<2; cath++) {
454 cnew.fX[cath]=Float_t(xm[ico][1]);
455 cnew.fY[cath]=Float_t(ym[ico][0]);
456 cnew.fZ[cath]=fZPlane;
458 cnew.fMultiplicity[cath]=c->fMultiplicity[cath];
459 for (i=0; i<fMul[cath]; i++) {
460 cnew.fIndexMap[i][cath]=c->fIndexMap[i][cath];
461 fSeg[cath]->SetPad(fIx[i][cath], fIy[i][cath]);
463 fprintf(stderr,"\nRawCluster %d cath %d\n",ico,cath);
464 fprintf(stderr,"mult_av %d\n",c->fMultiplicity[cath]);
465 FillCluster(&cnew,cath);
467 cnew.fClusterType=cnew.PhysicsContribution();
476 // ******* iacc = 2 *******
477 // Two combinations found between the 2 cathodes
479 // Was the same maximum taken twice
480 if ((accepted[0]&&accepted[1]) || (accepted[2]&&accepted[3])) {
481 fprintf(stderr,"\n Maximum taken twice !!!\n");
483 // Have a try !! with that
484 if (accepted[0]&&accepted[3]) {
496 fprintf(stderr,"\n cas (2) CombiDoubleMathiesonFit(c)\n");
497 chi2=CombiDoubleMathiesonFit(c);
498 // Int_t ndf = fgNbins[0]+fgNbins[1]-6;
499 // Float_t prob = TMath::Prob(chi2,ndf);
500 // prob2->Fill(prob);
501 // chi2_2->Fill(chi2);
505 // No ghosts ! No Problems ! - Perform one fit only !
506 if (accepted[0]&&accepted[3]) {
518 fprintf(stderr,"\n cas (2) CombiDoubleMathiesonFit(c)\n");
519 chi2=CombiDoubleMathiesonFit(c);
520 // Int_t ndf = fgNbins[0]+fgNbins[1]-6;
521 // Float_t prob = TMath::Prob(chi2,ndf);
522 // prob2->Fill(prob);
523 // chi2_2->Fill(chi2);
525 fprintf(stderr," chi2 %f\n",chi2);
529 // ******* iacc = 4 *******
530 // Four combinations found between the 2 cathodes
532 } else if (iacc==4) {
533 // Perform fits for the two possibilities !!
534 // Accept if charges are compatible on both cathodes
535 // If none are compatible, keep everything
541 fprintf(stderr,"\n cas (2) CombiDoubleMathiesonFit(c)\n");
542 chi2=CombiDoubleMathiesonFit(c);
543 // Int_t ndf = fgNbins[0]+fgNbins[1]-6;
544 // Float_t prob = TMath::Prob(chi2,ndf);
545 // prob2->Fill(prob);
546 // chi2_2->Fill(chi2);
548 fprintf(stderr," chi2 %f\n",chi2);
549 // store results of fit and postpone decision
550 Double_t sXFit[2],sYFit[2],sQrFit[2];
552 for (Int_t i=0;i<2;i++) {
563 fprintf(stderr,"\n cas (2) CombiDoubleMathiesonFit(c)\n");
564 chi2=CombiDoubleMathiesonFit(c);
565 // ndf = fgNbins[0]+fgNbins[1]-6;
566 // prob = TMath::Prob(chi2,ndf);
567 // prob2->Fill(prob);
568 // chi2_2->Fill(chi2);
570 fprintf(stderr," chi2 %f\n",chi2);
571 // We have all informations to perform the decision
572 // Compute the chi2 for the 2 possibilities
573 Float_t chi2fi,chi2si,chi2f,chi2s;
575 chi2f = (TMath::Log(fInput->TotalCharge(0)*fQrFit[0]
576 / (fInput->TotalCharge(1)*fQrFit[1]) )
577 / fInput->Response()->ChargeCorrel() );
579 chi2fi = (TMath::Log(fInput->TotalCharge(0)*(1-fQrFit[0])
580 / (fInput->TotalCharge(1)*(1-fQrFit[1])) )
581 / fInput->Response()->ChargeCorrel() );
582 chi2f += chi2fi*chi2fi;
584 chi2s = (TMath::Log(fInput->TotalCharge(0)*sQrFit[0]
585 / (fInput->TotalCharge(1)*sQrFit[1]) )
586 / fInput->Response()->ChargeCorrel() );
588 chi2si = (TMath::Log(fInput->TotalCharge(0)*(1-sQrFit[0])
589 / (fInput->TotalCharge(1)*(1-sQrFit[1])) )
590 / fInput->Response()->ChargeCorrel() );
591 chi2s += chi2si*chi2si;
593 // usefull to store the charge matching chi2 in the cluster
594 // fChi2[0]=sChi2[1]=chi2f;
595 // fChi2[1]=sChi2[0]=chi2s;
597 if (chi2f<=fGhostChi2Cut && chi2s<=fGhostChi2Cut)
599 if (chi2f>fGhostChi2Cut && chi2s>fGhostChi2Cut) {
605 if (chi2f<=fGhostChi2Cut)
607 if (chi2s<=fGhostChi2Cut) {
608 // retreive saved values
609 for (Int_t i=0;i<2;i++) {
620 } else if (fNLocal[0]==2 && fNLocal[1]==1) {
621 // ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
622 // (3) Two local maxima on cathode 1 and one maximum on cathode 2
623 // ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
625 Float_t xm[4][2], ym[4][2];
626 Float_t dpx, dpy, dx, dy;
627 Int_t ixm[4][2], iym[4][2];
628 Int_t isec, im1, ico;
630 // Form the 2x2 combinations
631 // 0-0, 0-1, 1-0, 1-1
633 for (im1=0; im1<2; im1++) {
634 xm[ico][0]=fX[fIndLocal[im1][0]][0];
635 ym[ico][0]=fY[fIndLocal[im1][0]][0];
636 xm[ico][1]=fX[fIndLocal[0][1]][1];
637 ym[ico][1]=fY[fIndLocal[0][1]][1];
639 ixm[ico][0]=fIx[fIndLocal[im1][0]][0];
640 iym[ico][0]=fIy[fIndLocal[im1][0]][0];
641 ixm[ico][1]=fIx[fIndLocal[0][1]][1];
642 iym[ico][1]=fIy[fIndLocal[0][1]][1];
645 // ico = 0 : first local maximum on cathodes 1 and 2
646 // ico = 1 : second local maximum on cathode 1 and first on cathode 2
648 // Analyse the combinations and keep those that are possible !
649 // For each combination check consistency in x and y
653 // In case of staggering maxima are displaced by exactly half the pad-size in y.
654 // We have to take into account the numerical precision in the consistency check;
658 for (ico=0; ico<2; ico++) {
659 accepted[ico]=kFALSE;
660 isec=fSeg[0]->Sector(ixm[ico][0], iym[ico][0]);
661 dpx=fSeg[0]->Dpx(isec)/2.;
662 dx=TMath::Abs(xm[ico][0]-xm[ico][1]);
663 isec=fSeg[1]->Sector(ixm[ico][1], iym[ico][1]);
664 dpy=fSeg[1]->Dpy(isec)/2.;
665 dy=TMath::Abs(ym[ico][0]-ym[ico][1]);
667 printf("\n %i %f %f %f %f \n", ico, ym[ico][0], ym[ico][1], dy, dpy );
668 if ((dx <= dpx) && (dy <= dpy+eps)) {
674 accepted[ico]=kFALSE;
682 // Initial value for charge ratios
683 fQrInit[0]=Float_t(fQ[fIndLocal[0][0]][0])/
684 Float_t(fQ[fIndLocal[0][0]][0]+fQ[fIndLocal[1][0]][0]);
685 fQrInit[1]=fQrInit[0];
687 if (accepted[0] && accepted[1]) {
689 fXInit[0]=0.5*(xm[0][1]+xm[0][0]);
691 fXInit[1]=0.5*(xm[0][1]+xm[1][0]);
695 chi23=CombiDoubleMathiesonFit(c);
704 } else if (accepted[0]) {
709 chi21=CombiDoubleMathiesonFit(c);
710 // Int_t ndf = fgNbins[0]+fgNbins[1]-6;
711 // Float_t prob = TMath::Prob(chi2,ndf);
712 // prob2->Fill(prob);
713 // chi2_2->Fill(chi21);
715 fprintf(stderr," chi2 %f\n",chi21);
716 if (chi21<10) Split(c);
717 } else if (accepted[1]) {
722 chi22=CombiDoubleMathiesonFit(c);
723 // Int_t ndf = fgNbins[0]+fgNbins[1]-6;
724 // Float_t prob = TMath::Prob(chi2,ndf);
725 // prob2->Fill(prob);
726 // chi2_2->Fill(chi22);
728 fprintf(stderr," chi2 %f\n",chi22);
729 if (chi22<10) Split(c);
732 if (chi21 > 10 && chi22 > 10 && chi23 > 10) {
733 // We keep only the combination found (X->cathode 2, Y->cathode 1)
734 for (Int_t ico=0; ico<2; ico++) {
736 AliMUONRawCluster cnew;
738 for (cath=0; cath<2; cath++) {
739 cnew.fX[cath]=Float_t(xm[ico][1]);
740 cnew.fY[cath]=Float_t(ym[ico][0]);
741 cnew.fZ[cath]=fZPlane;
742 cnew.fMultiplicity[cath]=c->fMultiplicity[cath];
743 for (i=0; i<fMul[cath]; i++) {
744 cnew.fIndexMap[i][cath]=c->fIndexMap[i][cath];
745 fSeg[cath]->SetPad(fIx[i][cath], fIy[i][cath]);
747 fprintf(stderr,"\nRawCluster %d cath %d\n",ico,cath);
748 fprintf(stderr,"mult_av %d\n",c->fMultiplicity[cath]);
749 FillCluster(&cnew,cath);
751 cnew.fClusterType=cnew.PhysicsContribution();
758 // +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
759 // (3') One local maximum on cathode 1 and two maxima on cathode 2
760 // +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
761 } else if (fNLocal[0]==1 && fNLocal[1]==2) {
762 Float_t xm[4][2], ym[4][2];
763 Float_t dpx, dpy, dx, dy;
764 Int_t ixm[4][2], iym[4][2];
765 Int_t isec, im1, ico;
767 // Form the 2x2 combinations
768 // 0-0, 0-1, 1-0, 1-1
770 for (im1=0; im1<2; im1++) {
771 xm[ico][0]=fX[fIndLocal[0][0]][0];
772 ym[ico][0]=fY[fIndLocal[0][0]][0];
773 xm[ico][1]=fX[fIndLocal[im1][1]][1];
774 ym[ico][1]=fY[fIndLocal[im1][1]][1];
776 ixm[ico][0]=fIx[fIndLocal[0][0]][0];
777 iym[ico][0]=fIy[fIndLocal[0][0]][0];
778 ixm[ico][1]=fIx[fIndLocal[im1][1]][1];
779 iym[ico][1]=fIy[fIndLocal[im1][1]][1];
782 // ico = 0 : first local maximum on cathodes 1 and 2
783 // ico = 1 : first local maximum on cathode 1 and second on cathode 2
785 // Analyse the combinations and keep those that are possible !
786 // For each combination check consistency in x and y
790 // In case of staggering maxima are displaced by exactly half the pad-size in y.
791 // We have to take into account the numerical precision in the consistency check;
795 for (ico=0; ico<2; ico++) {
796 accepted[ico]=kFALSE;
797 isec=fSeg[0]->Sector(ixm[ico][0], iym[ico][0]);
798 dpx=fSeg[0]->Dpx(isec)/2.;
799 dx=TMath::Abs(xm[ico][0]-xm[ico][1]);
800 isec=fSeg[1]->Sector(ixm[ico][1], iym[ico][1]);
801 dpy=fSeg[1]->Dpy(isec)/2.;
802 dy=TMath::Abs(ym[ico][0]-ym[ico][1]);
804 printf("\n %i %f %f %f %f \n", ico, ym[ico][0], ym[ico][1], dy, dpy );
805 if ((dx <= dpx) && (dy <= dpy+eps)) {
808 fprintf(stderr,"ico %d\n",ico);
812 accepted[ico]=kFALSE;
820 fQrInit[1]=Float_t(fQ[fIndLocal[0][1]][1])/
821 Float_t(fQ[fIndLocal[0][1]][1]+fQ[fIndLocal[1][1]][1]);
823 fQrInit[0]=fQrInit[1];
826 if (accepted[0] && accepted[1]) {
828 fYInit[0]=0.5*(ym[0][0]+ym[0][1]);
830 fYInit[1]=0.5*(ym[0][0]+ym[1][1]);
833 chi23=CombiDoubleMathiesonFit(c);
842 } else if (accepted[0]) {
847 chi21=CombiDoubleMathiesonFit(c);
848 // Int_t ndf = fgNbins[0]+fgNbins[1]-6;
849 // Float_t prob = TMath::Prob(chi2,ndf);
850 // prob2->Fill(prob);
851 // chi2_2->Fill(chi21);
853 fprintf(stderr," chi2 %f\n",chi21);
854 if (chi21<10) Split(c);
855 } else if (accepted[1]) {
860 chi22=CombiDoubleMathiesonFit(c);
861 // Int_t ndf = fgNbins[0]+fgNbins[1]-6;
862 // Float_t prob = TMath::Prob(chi2,ndf);
863 // prob2->Fill(prob);
864 // chi2_2->Fill(chi22);
866 fprintf(stderr," chi2 %f\n",chi22);
867 if (chi22<10) Split(c);
870 if (chi21 > 10 && chi22 > 10 && chi23 > 10) {
871 //We keep only the combination found (X->cathode 2, Y->cathode 1)
872 for (Int_t ico=0; ico<2; ico++) {
874 AliMUONRawCluster cnew;
876 for (cath=0; cath<2; cath++) {
877 cnew.fX[cath]=Float_t(xm[ico][1]);
878 cnew.fY[cath]=Float_t(ym[ico][0]);
879 cnew.fZ[cath]=fZPlane;
880 cnew.fMultiplicity[cath]=c->fMultiplicity[cath];
881 for (i=0; i<fMul[cath]; i++) {
882 cnew.fIndexMap[i][cath]=c->fIndexMap[i][cath];
883 fSeg[cath]->SetPad(fIx[i][cath], fIy[i][cath]);
885 fprintf(stderr,"\nRawCluster %d cath %d\n",ico,cath);
886 fprintf(stderr,"mult_av %d\n",c->fMultiplicity[cath]);
887 FillCluster(&cnew,cath);
889 cnew.fClusterType=cnew.PhysicsContribution();
896 // ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
897 // (4) At least three local maxima on cathode 1 or on cathode 2
898 // ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
899 } else if (fNLocal[0]>2 || fNLocal[1]>2) {
900 Int_t param = fNLocal[0]*fNLocal[1];
903 Float_t ** xm = new Float_t * [param];
904 for (ii=0; ii<param; ii++) xm[ii]=new Float_t [2];
905 Float_t ** ym = new Float_t * [param];
906 for (ii=0; ii<param; ii++) ym[ii]=new Float_t [2];
907 Int_t ** ixm = new Int_t * [param];
908 for (ii=0; ii<param; ii++) ixm[ii]=new Int_t [2];
909 Int_t ** iym = new Int_t * [param];
910 for (ii=0; ii<param; ii++) iym[ii]=new Int_t [2];
913 Float_t dpx, dpy, dx, dy;
916 for (Int_t im1=0; im1<fNLocal[0]; im1++) {
917 for (Int_t im2=0; im2<fNLocal[1]; im2++) {
918 xm[ico][0]=fX[fIndLocal[im1][0]][0];
919 ym[ico][0]=fY[fIndLocal[im1][0]][0];
920 xm[ico][1]=fX[fIndLocal[im2][1]][1];
921 ym[ico][1]=fY[fIndLocal[im2][1]][1];
923 ixm[ico][0]=fIx[fIndLocal[im1][0]][0];
924 iym[ico][0]=fIy[fIndLocal[im1][0]][0];
925 ixm[ico][1]=fIx[fIndLocal[im2][1]][1];
926 iym[ico][1]=fIy[fIndLocal[im2][1]][1];
933 fprintf(stderr,"nIco %d\n",nIco);
934 for (ico=0; ico<nIco; ico++) {
936 fprintf(stderr,"ico = %d\n",ico);
937 isec=fSeg[0]->Sector(ixm[ico][0], iym[ico][0]);
938 dpx=fSeg[0]->Dpx(isec)/2.;
939 dx=TMath::Abs(xm[ico][0]-xm[ico][1]);
940 isec=fSeg[1]->Sector(ixm[ico][1], iym[ico][1]);
941 dpy=fSeg[1]->Dpy(isec)/2.;
942 dy=TMath::Abs(ym[ico][0]-ym[ico][1]);
944 fprintf(stderr,"dx %f dpx %f dy %f dpy %f\n",dx,dpx,dy,dpy);
945 fprintf(stderr," X %f Y %f\n",xm[ico][1],ym[ico][0]);
947 if ((dx <= dpx) && (dy <= dpy)) {
949 fprintf(stderr,"ok\n");
951 AliMUONRawCluster cnew;
952 for (cath=0; cath<2; cath++) {
953 cnew.fX[cath]=Float_t(xm[ico][1]);
954 cnew.fY[cath]=Float_t(ym[ico][0]);
955 cnew.fZ[cath]=fZPlane;
956 cnew.fMultiplicity[cath]=c->fMultiplicity[cath];
957 for (i=0; i<fMul[cath]; i++) {
958 cnew.fIndexMap[i][cath]=c->fIndexMap[i][cath];
959 fSeg[cath]->SetPad(fIx[i][cath], fIy[i][cath]);
961 FillCluster(&cnew,cath);
963 cnew.fClusterType=cnew.PhysicsContribution();
975 void AliMUONClusterFinderVS::FindLocalMaxima(AliMUONRawCluster* /*c*/)
977 // Find all local maxima of a cluster
979 printf("\n Find Local maxima !");
983 Int_t cath, cath1; // loops over cathodes
984 Int_t i; // loops over digits
985 Int_t j; // loops over cathodes
989 // counters for number of local maxima
990 fNLocal[0]=fNLocal[1]=0;
991 // flags digits as local maximum
992 Bool_t isLocal[100][2];
993 for (i=0; i<100;i++) {
994 isLocal[i][0]=isLocal[i][1]=kFALSE;
996 // number of next neighbours and arrays to store them
999 // loop over cathodes
1000 for (cath=0; cath<2; cath++) {
1001 // loop over cluster digits
1002 for (i=0; i<fMul[cath]; i++) {
1003 // get neighbours for that digit and assume that it is local maximum
1004 fSeg[cath]->Neighbours(fIx[i][cath], fIy[i][cath], &nn, x, y);
1005 isLocal[i][cath]=kTRUE;
1006 Int_t isec= fSeg[cath]->Sector(fIx[i][cath], fIy[i][cath]);
1007 Float_t a0 = fSeg[cath]->Dpx(isec)*fSeg[cath]->Dpy(isec);
1008 // loop over next neighbours, if at least one neighbour has higher charger assumption
1009 // digit is not local maximum
1010 for (j=0; j<nn; j++) {
1011 if (fHitMap[cath]->TestHit(x[j], y[j])==kEmpty) continue;
1012 digt=(AliMUONDigit*) fHitMap[cath]->GetHit(x[j], y[j]);
1013 isec=fSeg[cath]->Sector(x[j], y[j]);
1014 Float_t a1 = fSeg[cath]->Dpx(isec)*fSeg[cath]->Dpy(isec);
1015 if (digt->Signal()/a1 > fQ[i][cath]/a0) {
1016 isLocal[i][cath]=kFALSE;
1019 // handle special case of neighbouring pads with equal signal
1020 } else if (digt->Signal() == fQ[i][cath]) {
1021 if (fNLocal[cath]>0) {
1022 for (Int_t k=0; k<fNLocal[cath]; k++) {
1023 if (x[j]==fIx[fIndLocal[k][cath]][cath]
1024 && y[j]==fIy[fIndLocal[k][cath]][cath])
1026 isLocal[i][cath]=kFALSE;
1028 } // loop over local maxima
1029 } // are there already local maxima
1031 } // loop over next neighbours
1032 if (isLocal[i][cath]) {
1033 fIndLocal[fNLocal[cath]][cath]=i;
1036 } // loop over all digits
1037 } // loop over cathodes
1040 printf("\n Found %d %d %d %d local Maxima\n",
1041 fNLocal[0], fNLocal[1], fMul[0], fMul[1]);
1042 fprintf(stderr,"\n Cathode 1 local Maxima %d Multiplicite %d\n",fNLocal[0], fMul[0]);
1043 fprintf(stderr," Cathode 2 local Maxima %d Multiplicite %d\n",fNLocal[1], fMul[1]);
1049 if (fNLocal[1]==2 && (fNLocal[0]==1 || fNLocal[0]==0)) {
1050 Int_t iback=fNLocal[0];
1052 // Two local maxima on cathode 2 and one maximum on cathode 1
1053 // Look for local maxima considering up and down neighbours on the 1st cathode only
1055 // Loop over cluster digits
1059 for (i=0; i<fMul[cath]; i++) {
1060 isec=fSeg[cath]->Sector(fIx[i][cath],fIy[i][cath]);
1061 dpy=fSeg[cath]->Dpy(isec);
1062 dpx=fSeg[cath]->Dpx(isec);
1063 if (isLocal[i][cath]) continue;
1064 // Pad position should be consistent with position of local maxima on the opposite cathode
1065 if ((TMath::Abs(fX[i][cath]-fX[fIndLocal[0][cath1]][cath1]) > dpx/2.) &&
1066 (TMath::Abs(fX[i][cath]-fX[fIndLocal[1][cath1]][cath1]) > dpx/2.))
1069 // get neighbours for that digit and assume that it is local maximum
1070 isLocal[i][cath]=kTRUE;
1071 // compare signal to that on the two neighbours on the left and on the right
1072 // iNN counts the number of neighbours with signal, it should be 1 or 2
1076 ->FirstPad(fX[i][cath], fY[i][cath], fZPlane, 0., dpy);
1082 ix = fSeg[cath]->Ix();
1083 iy = fSeg[cath]->Iy();
1084 // skip the current pad
1085 if (iy == fIy[i][cath]) continue;
1087 if (fHitMap[cath]->TestHit(ix, iy)!=kEmpty) {
1089 digt=(AliMUONDigit*) fHitMap[cath]->GetHit(ix,iy);
1090 if (digt->Signal() > fQ[i][cath]) isLocal[i][cath]=kFALSE;
1092 } // Loop over pad neighbours in y
1093 if (isLocal[i][cath] && iNN>0) {
1094 fIndLocal[fNLocal[cath]][cath]=i;
1097 } // loop over all digits
1098 // if one additional maximum has been found we are happy
1099 // if more maxima have been found restore the previous situation
1102 "\n New search gives %d local maxima for cathode 1 \n",
1105 " %d local maxima for cathode 2 \n",
1108 if (fNLocal[cath]>2) {
1109 fNLocal[cath]=iback;
1112 } // 1,2 local maxima
1114 if (fNLocal[0]==2 && (fNLocal[1]==1 || fNLocal[1]==0)) {
1115 Int_t iback=fNLocal[1];
1117 // Two local maxima on cathode 1 and one maximum on cathode 2
1118 // Look for local maxima considering left and right neighbours on the 2nd cathode only
1121 Float_t eps = 1.e-5;
1124 // Loop over cluster digits
1125 for (i=0; i<fMul[cath]; i++) {
1126 isec=fSeg[cath]->Sector(fIx[i][cath],fIy[i][cath]);
1127 dpx=fSeg[cath]->Dpx(isec);
1128 dpy=fSeg[cath]->Dpy(isec);
1129 if (isLocal[i][cath]) continue;
1130 // Pad position should be consistent with position of local maxima on the opposite cathode
1131 if ((TMath::Abs(fY[i][cath]-fY[fIndLocal[0][cath1]][cath1]) > dpy/2.+eps) &&
1132 (TMath::Abs(fY[i][cath]-fY[fIndLocal[1][cath1]][cath1]) > dpy/2.+eps))
1136 // get neighbours for that digit and assume that it is local maximum
1137 isLocal[i][cath]=kTRUE;
1138 // compare signal to that on the two neighbours on the left and on the right
1140 // iNN counts the number of neighbours with signal, it should be 1 or 2
1143 ->FirstPad(fX[i][cath], fY[i][cath], fZPlane, dpx, 0.);
1150 ix = fSeg[cath]->Ix();
1151 iy = fSeg[cath]->Iy();
1153 // skip the current pad
1154 if (ix == fIx[i][cath]) continue;
1156 if (fHitMap[cath]->TestHit(ix, iy)!=kEmpty) {
1158 digt=(AliMUONDigit*) fHitMap[cath]->GetHit(ix,iy);
1159 if (digt->Signal() > fQ[i][cath]) isLocal[i][cath]=kFALSE;
1161 } // Loop over pad neighbours in x
1162 if (isLocal[i][cath] && iNN>0) {
1163 fIndLocal[fNLocal[cath]][cath]=i;
1166 } // loop over all digits
1167 // if one additional maximum has been found we are happy
1168 // if more maxima have been found restore the previous situation
1170 fprintf(stderr,"\n New search gives %d local maxima for cathode 1 \n",fNLocal[0]);
1171 fprintf(stderr,"\n %d local maxima for cathode 2 \n",fNLocal[1]);
1172 printf("\n New search gives %d %d \n",fNLocal[0],fNLocal[1]);
1174 if (fNLocal[cath]>2) {
1175 fNLocal[cath]=iback;
1177 } // 2,1 local maxima
1181 void AliMUONClusterFinderVS::FillCluster(AliMUONRawCluster* c, Int_t flag, Int_t cath)
1184 // Completes cluster information starting from list of digits
1191 c->fPeakSignal[cath]=c->fPeakSignal[0];
1193 c->fPeakSignal[cath]=0;
1204 fprintf(stderr,"\n fPeakSignal %d\n",c->fPeakSignal[cath]);
1205 for (Int_t i=0; i<c->fMultiplicity[cath]; i++)
1207 dig= fInput->Digit(cath,c->fIndexMap[i][cath]);
1208 ix=dig->PadX()+c->fOffsetMap[i][cath];
1210 Int_t q=dig->Signal();
1211 if (!flag) q=Int_t(q*c->fContMap[i][cath]);
1212 // fprintf(stderr,"q %d c->fPeakSignal[ %d ] %d\n",q,cath,c->fPeakSignal[cath]);
1213 if (dig->Physics() >= dig->Signal()) {
1214 c->fPhysicsMap[i]=2;
1215 } else if (dig->Physics() == 0) {
1216 c->fPhysicsMap[i]=0;
1217 } else c->fPhysicsMap[i]=1;
1221 fprintf(stderr,"q %d c->fPeakSignal[cath] %d\n",q,c->fPeakSignal[cath]);
1222 // peak signal and track list
1223 if (q>c->fPeakSignal[cath]) {
1224 c->fPeakSignal[cath]=q;
1225 c->fTracks[0]=dig->Hit();
1226 c->fTracks[1]=dig->Track(0);
1227 c->fTracks[2]=dig->Track(1);
1228 // fprintf(stderr," c->fTracks[0] %d c->fTracks[1] %d\n",dig->fHit,dig->fTracks[0]);
1232 fSeg[cath]->GetPadC(ix, iy, x, y, z);
1237 } // loop over digits
1239 fprintf(stderr," fin du cluster c\n");
1243 c->fX[cath]/=c->fQ[cath];
1245 c->fX[cath]=fSeg[cath]->GetAnod(c->fX[cath]);
1246 c->fY[cath]/=c->fQ[cath];
1248 // apply correction to the coordinate along the anode wire
1252 fSeg[cath]->GetPadI(x, y, fZPlane, ix, iy);
1253 fSeg[cath]->GetPadC(ix, iy, x, y, z);
1254 Int_t isec=fSeg[cath]->Sector(ix,iy);
1255 TF1* cogCorr = fSeg[cath]->CorrFunc(isec-1);
1258 Float_t yOnPad=(c->fY[cath]-y)/fSeg[cath]->Dpy(isec);
1259 c->fY[cath]=c->fY[cath]-cogCorr->Eval(yOnPad, 0, 0);
1264 void AliMUONClusterFinderVS::FillCluster(AliMUONRawCluster* c, Int_t cath)
1267 // Completes cluster information starting from list of digits
1277 Float_t xpad, ypad, zpad;
1280 for (Int_t i=0; i<c->fMultiplicity[cath]; i++)
1282 dig = fInput->Digit(cath,c->fIndexMap[i][cath]);
1284 GetPadC(dig->PadX(),dig->PadY(),xpad,ypad, zpad);
1286 fprintf(stderr,"x %f y %f cx %f cy %f\n",xpad,ypad,c->fX[0],c->fY[0]);
1287 dx = xpad - c->fX[0];
1288 dy = ypad - c->fY[0];
1289 dr = TMath::Sqrt(dx*dx+dy*dy);
1294 fprintf(stderr," dr %f\n",dr);
1295 Int_t q=dig->Signal();
1296 if (dig->Physics() >= dig->Signal()) {
1297 c->fPhysicsMap[i]=2;
1298 } else if (dig->Physics() == 0) {
1299 c->fPhysicsMap[i]=0;
1300 } else c->fPhysicsMap[i]=1;
1301 c->fPeakSignal[cath]=q;
1302 c->fTracks[0]=dig->Hit();
1303 c->fTracks[1]=dig->Track(0);
1304 c->fTracks[2]=dig->Track(1);
1306 fprintf(stderr," c->fTracks[0] %d c->fTracks[1] %d\n",dig->Hit(),
1310 } // loop over digits
1312 // apply correction to the coordinate along the anode wire
1314 c->fX[cath]=fSeg[cath]->GetAnod(c->fX[cath]);
1317 void AliMUONClusterFinderVS::FindCluster(Int_t i, Int_t j, Int_t cath, AliMUONRawCluster &c){
1321 // Find a super cluster on both cathodes
1324 // Add i,j as element of the cluster
1327 Int_t idx = fHitMap[cath]->GetHitIndex(i,j);
1328 AliMUONDigit* dig = (AliMUONDigit*) fHitMap[cath]->GetHit(i,j);
1329 Int_t q=dig->Signal();
1330 Int_t theX=dig->PadX();
1331 Int_t theY=dig->PadY();
1333 if (q > TMath::Abs(c.fPeakSignal[0]) && q > TMath::Abs(c.fPeakSignal[1])) {
1334 c.fPeakSignal[cath]=q;
1335 c.fTracks[0]=dig->Hit();
1336 c.fTracks[1]=dig->Track(0);
1337 c.fTracks[2]=dig->Track(1);
1341 // Make sure that list of digits is ordered
1343 Int_t mu=c.fMultiplicity[cath];
1344 c.fIndexMap[mu][cath]=idx;
1346 if (dig->Physics() >= dig->Signal()) {
1347 c.fPhysicsMap[mu]=2;
1348 } else if (dig->Physics() == 0) {
1349 c.fPhysicsMap[mu]=0;
1350 } else c.fPhysicsMap[mu]=1;
1354 for (Int_t ind = mu-1; ind >= 0; ind--) {
1355 Int_t ist=(c.fIndexMap)[ind][cath];
1356 Int_t ql=fInput->Digit(cath, ist)->Signal();
1357 Int_t ix=fInput->Digit(cath, ist)->PadX();
1358 Int_t iy=fInput->Digit(cath, ist)->PadY();
1360 if (q>ql || (q==ql && theX > ix && theY < iy)) {
1361 c.fIndexMap[ind][cath]=idx;
1362 c.fIndexMap[ind+1][cath]=ist;
1370 c.fMultiplicity[cath]++;
1371 if (c.fMultiplicity[cath] >= 50 ) {
1372 printf("FindCluster - multiplicity >50 %d \n",c.fMultiplicity[0]);
1373 c.fMultiplicity[cath]=49;
1376 // Prepare center of gravity calculation
1378 fSeg[cath]->GetPadC(i, j, x, y, z);
1384 // Flag hit as "taken"
1385 fHitMap[cath]->FlagHit(i,j);
1387 // Now look recursively for all neighbours and pad hit on opposite cathode
1389 // Loop over neighbours
1393 Int_t xList[10], yList[10];
1394 fSeg[cath]->Neighbours(i,j,&nn,xList,yList);
1395 for (Int_t in=0; in<nn; in++) {
1399 if (fHitMap[cath]->TestHit(ix,iy)==kUnused) {
1401 printf("\n Neighbours %d %d %d", cath, ix, iy);
1402 FindCluster(ix, iy, cath, c);
1407 Int_t iXopp[50], iYopp[50];
1409 // Neighbours on opposite cathode
1410 // Take into account that several pads can overlap with the present pad
1411 Int_t isec=fSeg[cath]->Sector(i,j);
1417 dx = (fSeg[cath]->Dpx(isec))/2.;
1422 dy = (fSeg[cath]->Dpy(isec))/2;
1424 // loop over pad neighbours on opposite cathode
1425 for (fSeg[iop]->FirstPad(x, y, fZPlane, dx, dy);
1426 fSeg[iop]->MorePads();
1427 fSeg[iop]->NextPad())
1430 ix = fSeg[iop]->Ix(); iy = fSeg[iop]->Iy();
1431 if (fDebugLevel > 1)
1432 printf("\n ix, iy: %f %f %f %d %d %d", x,y,z,ix, iy, fSector);
1433 if (fHitMap[iop]->TestHit(ix,iy)==kUnused){
1436 if (fDebugLevel > 1)
1437 printf("\n Opposite %d %d %d", iop, ix, iy);
1440 } // Loop over pad neighbours
1441 // This had to go outside the loop since recursive calls inside the iterator are not possible
1444 for (jopp=0; jopp<nOpp; jopp++) {
1445 if (fHitMap[iop]->TestHit(iXopp[jopp],iYopp[jopp]) == kUnused)
1446 FindCluster(iXopp[jopp], iYopp[jopp], iop, c);
1450 //_____________________________________________________________________________
1452 void AliMUONClusterFinderVS::FindRawClusters()
1455 // MUON cluster finder from digits -- finds neighbours on both cathodes and
1456 // fills the tree with raw clusters
1460 // Return if no input datad available
1461 if (!fInput->NDigits(0) && !fInput->NDigits(1)) return;
1463 fSeg[0] = fInput->Segmentation(0);
1464 fSeg[1] = fInput->Segmentation(1);
1466 fHitMap[0] = new AliMUONHitMapA1(fSeg[0], fInput->Digits(0));
1467 fHitMap[1] = new AliMUONHitMapA1(fSeg[1], fInput->Digits(1));
1475 fHitMap[0]->FillHits();
1476 fHitMap[1]->FillHits();
1478 // Outer Loop over Cathodes
1479 for (cath=0; cath<2; cath++) {
1480 for (ndig=0; ndig<fInput->NDigits(cath); ndig++) {
1481 dig = fInput->Digit(cath, ndig);
1482 Int_t i=dig->PadX();
1483 Int_t j=dig->PadY();
1484 if (fHitMap[cath]->TestHit(i,j)==kUsed ||fHitMap[0]->TestHit(i,j)==kEmpty) {
1489 fprintf(stderr,"\n CATHODE %d CLUSTER %d\n",cath,ncls);
1490 AliMUONRawCluster c;
1491 c.fMultiplicity[0]=0;
1492 c.fMultiplicity[1]=0;
1493 c.fPeakSignal[cath]=dig->Signal();
1494 c.fTracks[0]=dig->Hit();
1495 c.fTracks[1]=dig->Track(0);
1496 c.fTracks[2]=dig->Track(1);
1497 // tag the beginning of cluster list in a raw cluster
1500 fSeg[cath]->GetPadC(i,j,xcu, ycu, fZPlane);
1501 fSector= fSeg[cath]->Sector(i,j)/100;
1503 printf("\n New Seed %d %d ", i,j);
1506 FindCluster(i,j,cath,c);
1507 // ^^^^^^^^^^^^^^^^^^^^^^^^
1508 // center of gravity
1509 if (c.fX[0]!=0.) c.fX[0] /= c.fQ[0];
1511 c.fX[0]=fSeg[0]->GetAnod(c.fX[0]);
1512 if (c.fY[0]!=0.) c.fY[0] /= c.fQ[0];
1514 if(c.fQ[1]!=0.) c.fX[1] /= c.fQ[1];
1517 c.fX[1]=fSeg[0]->GetAnod(c.fX[1]);
1518 if(c.fQ[1]!=0.) c.fY[1] /= c.fQ[1];
1524 fprintf(stderr,"\n Cathode 1 multiplicite %d X(CG) %f Y(CG) %f\n",
1525 c.fMultiplicity[0],c.fX[0],c.fY[0]);
1526 fprintf(stderr," Cathode 2 multiplicite %d X(CG) %f Y(CG) %f\n",
1527 c.fMultiplicity[1],c.fX[1],c.fY[1]);
1529 // Analyse cluster and decluster if necessary
1532 c.fNcluster[1]=fNRawClusters;
1533 c.fClusterType=c.PhysicsContribution();
1540 // reset Cluster object
1541 { // begin local scope
1542 for (int k=0;k<c.fMultiplicity[0];k++) c.fIndexMap[k][0]=0;
1543 } // end local scope
1545 { // begin local scope
1546 for (int k=0;k<c.fMultiplicity[1];k++) c.fIndexMap[k][1]=0;
1547 } // end local scope
1549 c.fMultiplicity[0]=c.fMultiplicity[0]=0;
1553 } // end loop cathodes
1558 Float_t AliMUONClusterFinderVS::SingleMathiesonFit(AliMUONRawCluster *c, Int_t cath)
1560 // Performs a single Mathieson fit on one cathode
1562 Double_t arglist[20];
1564 AliMUONClusterInput& clusterInput = *(AliMUONClusterInput::Instance());
1566 clusterInput.Fitter()->SetFCN(fcnS1);
1567 clusterInput.Fitter()->mninit(2,10,7);
1568 clusterInput.Fitter()->SetPrintLevel(-1+fDebugLevel);
1570 clusterInput.Fitter()->mnexcm("SET NOW", arglist, 0, ierflag);
1571 // Set starting values
1572 static Double_t vstart[2];
1577 // lower and upper limits
1578 static Double_t lower[2], upper[2];
1580 fSeg[cath]->GetPadI(c->fX[cath], c->fY[cath], fZPlane, ix, iy);
1581 Int_t isec=fSeg[cath]->Sector(ix, iy);
1582 lower[0]=vstart[0]-fSeg[cath]->Dpx(isec)/2;
1583 lower[1]=vstart[1]-fSeg[cath]->Dpy(isec)/2;
1585 upper[0]=lower[0]+fSeg[cath]->Dpx(isec);
1586 upper[1]=lower[1]+fSeg[cath]->Dpy(isec);
1589 static Double_t step[2]={0.0005, 0.0005};
1591 clusterInput.Fitter()->mnparm(0,"x1",vstart[0],step[0],lower[0],upper[0],ierflag);
1592 clusterInput.Fitter()->mnparm(1,"y1",vstart[1],step[1],lower[1],upper[1],ierflag);
1593 // ready for minimisation
1597 clusterInput.Fitter()->mnexcm("SET NOGR", arglist, 0, ierflag);
1598 clusterInput.Fitter()->mnexcm("MIGRAD", arglist, 0, ierflag);
1599 // clusterInput.Fitter()->mnexcm("EXIT" , arglist, 0, ierflag);
1600 Double_t fmin, fedm, errdef;
1601 Int_t npari, nparx, istat;
1603 clusterInput.Fitter()->mnstat(fmin, fedm, errdef, npari, nparx, istat);
1607 // Get fitted parameters
1608 Double_t xrec, yrec;
1610 Double_t epxz, b1, b2;
1612 clusterInput.Fitter()->mnpout(0, chname, xrec, epxz, b1, b2, ierflg);
1613 clusterInput.Fitter()->mnpout(1, chname, yrec, epxz, b1, b2, ierflg);
1619 Float_t AliMUONClusterFinderVS::CombiSingleMathiesonFit(AliMUONRawCluster * /*c*/)
1621 // Perform combined Mathieson fit on both cathode planes
1623 Double_t arglist[20];
1625 AliMUONClusterInput& clusterInput = *(AliMUONClusterInput::Instance());
1626 clusterInput.Fitter()->SetFCN(fcnCombiS1);
1627 clusterInput.Fitter()->mninit(2,10,7);
1628 clusterInput.Fitter()->SetPrintLevel(-1+fDebugLevel);
1630 clusterInput.Fitter()->mnexcm("SET NOW", arglist, 0, ierflag);
1631 static Double_t vstart[2];
1632 vstart[0]=fXInit[0];
1633 vstart[1]=fYInit[0];
1636 // lower and upper limits
1637 static Float_t lower[2], upper[2];
1639 fSeg[0]->GetPadI(fXInit[0], fYInit[0], fZPlane, ix, iy);
1640 isec=fSeg[0]->Sector(ix, iy);
1641 Float_t dpy=fSeg[0]->Dpy(isec);
1642 fSeg[1]->GetPadI(fXInit[0], fYInit[0], fZPlane, ix, iy);
1643 isec=fSeg[1]->Sector(ix, iy);
1644 Float_t dpx=fSeg[1]->Dpx(isec);
1647 Float_t xdum, ydum, zdum;
1649 // Find save upper and lower limits
1653 for (fSeg[1]->FirstPad(fXInit[0], fYInit[0], fZPlane, dpx, 0.);
1654 fSeg[1]->MorePads(); fSeg[1]->NextPad())
1656 ix=fSeg[1]->Ix(); iy=fSeg[1]->Iy();
1657 fSeg[1]->GetPadC(ix,iy, upper[0], ydum, zdum);
1658 if (icount ==0) lower[0]=upper[0];
1662 if (lower[0]>upper[0]) {xdum=lower[0]; lower[0]=upper[0]; upper[0]=xdum;}
1666 printf("\n single y %f %f", fXInit[0], fYInit[0]);
1668 for (fSeg[0]->FirstPad(fXInit[0], fYInit[0], fZPlane, 0., dpy);
1669 fSeg[0]->MorePads(); fSeg[0]->NextPad())
1671 ix=fSeg[0]->Ix(); iy=fSeg[0]->Iy();
1672 fSeg[0]->GetPadC(ix,iy,xdum,upper[1],zdum);
1673 if (icount ==0) lower[1]=upper[1];
1676 printf("\n upper lower %d %f %f", icount, upper[1], lower[1]);
1679 if (lower[1]>upper[1]) {xdum=lower[1]; lower[1]=upper[1]; upper[1]=xdum;}
1682 static Double_t step[2]={0.00001, 0.0001};
1684 clusterInput.Fitter()->mnparm(0,"x1",vstart[0],step[0],lower[0],upper[0],ierflag);
1685 clusterInput.Fitter()->mnparm(1,"y1",vstart[1],step[1],lower[1],upper[1],ierflag);
1686 // ready for minimisation
1690 clusterInput.Fitter()->mnexcm("SET NOGR", arglist, 0, ierflag);
1691 clusterInput.Fitter()->mnexcm("MIGRAD", arglist, 0, ierflag);
1692 // clusterInput.Fitter()->mnexcm("EXIT" , arglist, 0, ierflag);
1693 Double_t fmin, fedm, errdef;
1694 Int_t npari, nparx, istat;
1696 clusterInput.Fitter()->mnstat(fmin, fedm, errdef, npari, nparx, istat);
1700 // Get fitted parameters
1701 Double_t xrec, yrec;
1703 Double_t epxz, b1, b2;
1705 clusterInput.Fitter()->mnpout(0, chname, xrec, epxz, b1, b2, ierflg);
1706 clusterInput.Fitter()->mnpout(1, chname, yrec, epxz, b1, b2, ierflg);
1712 Bool_t AliMUONClusterFinderVS::DoubleMathiesonFit(AliMUONRawCluster * /*c*/, Int_t cath)
1714 // Performs a double Mathieson fit on one cathode
1718 // Initialise global variables for fit
1719 Double_t arglist[20];
1721 AliMUONClusterInput& clusterInput = *(AliMUONClusterInput::Instance());
1722 clusterInput.Fitter()->SetFCN(fcnS2);
1723 clusterInput.Fitter()->mninit(5,10,7);
1724 clusterInput.Fitter()->SetPrintLevel(-1+fDebugLevel);
1726 clusterInput.Fitter()->mnexcm("SET NOW", arglist, 0, ierflag);
1727 // Set starting values
1728 static Double_t vstart[5];
1729 vstart[0]=fX[fIndLocal[0][cath]][cath];
1730 vstart[1]=fY[fIndLocal[0][cath]][cath];
1731 vstart[2]=fX[fIndLocal[1][cath]][cath];
1732 vstart[3]=fY[fIndLocal[1][cath]][cath];
1733 vstart[4]=Float_t(fQ[fIndLocal[0][cath]][cath])/
1734 Float_t(fQ[fIndLocal[0][cath]][cath]+fQ[fIndLocal[1][cath]][cath]);
1735 // lower and upper limits
1736 static Float_t lower[5], upper[5];
1737 Int_t isec=fSeg[cath]->Sector(fIx[fIndLocal[0][cath]][cath], fIy[fIndLocal[0][cath]][cath]);
1738 lower[0]=vstart[0]-fSeg[cath]->Dpx(isec);
1739 lower[1]=vstart[1]-fSeg[cath]->Dpy(isec);
1741 upper[0]=lower[0]+2.*fSeg[cath]->Dpx(isec);
1742 upper[1]=lower[1]+2.*fSeg[cath]->Dpy(isec);
1744 isec=fSeg[cath]->Sector(fIx[fIndLocal[1][cath]][cath], fIy[fIndLocal[1][cath]][cath]);
1745 lower[2]=vstart[2]-fSeg[cath]->Dpx(isec)/2;
1746 lower[3]=vstart[3]-fSeg[cath]->Dpy(isec)/2;
1748 upper[2]=lower[2]+fSeg[cath]->Dpx(isec);
1749 upper[3]=lower[3]+fSeg[cath]->Dpy(isec);
1754 static Double_t step[5]={0.0005, 0.0005, 0.0005, 0.0005, 0.0001};
1756 clusterInput.Fitter()->mnparm(0,"x1",vstart[0],step[0],lower[0],upper[0],ierflag);
1757 clusterInput.Fitter()->mnparm(1,"y1",vstart[1],step[1],lower[1],upper[1],ierflag);
1758 clusterInput.Fitter()->mnparm(2,"x2",vstart[2],step[2],lower[2],upper[2],ierflag);
1759 clusterInput.Fitter()->mnparm(3,"y2",vstart[3],step[3],lower[3],upper[3],ierflag);
1760 clusterInput.Fitter()->mnparm(4,"a0",vstart[4],step[4],lower[4],upper[4],ierflag);
1761 // ready for minimisation
1765 clusterInput.Fitter()->mnexcm("SET NOGR", arglist, 0, ierflag);
1766 clusterInput.Fitter()->mnexcm("MIGRAD", arglist, 0, ierflag);
1767 // clusterInput.Fitter()->mnexcm("EXIT" , arglist, 0, ierflag);
1768 // Get fitted parameters
1769 Double_t xrec[2], yrec[2], qfrac;
1771 Double_t epxz, b1, b2;
1773 clusterInput.Fitter()->mnpout(0, chname, xrec[0], epxz, b1, b2, ierflg);
1774 clusterInput.Fitter()->mnpout(1, chname, yrec[0], epxz, b1, b2, ierflg);
1775 clusterInput.Fitter()->mnpout(2, chname, xrec[1], epxz, b1, b2, ierflg);
1776 clusterInput.Fitter()->mnpout(3, chname, yrec[1], epxz, b1, b2, ierflg);
1777 clusterInput.Fitter()->mnpout(4, chname, qfrac, epxz, b1, b2, ierflg);
1779 Double_t fmin, fedm, errdef;
1780 Int_t npari, nparx, istat;
1782 clusterInput.Fitter()->mnstat(fmin, fedm, errdef, npari, nparx, istat);
1787 Float_t AliMUONClusterFinderVS::CombiDoubleMathiesonFit(AliMUONRawCluster * /*c*/)
1790 // Perform combined double Mathieson fit on both cathode planes
1792 Double_t arglist[20];
1794 AliMUONClusterInput& clusterInput = *(AliMUONClusterInput::Instance());
1795 clusterInput.Fitter()->SetFCN(fcnCombiS2);
1796 clusterInput.Fitter()->mninit(6,10,7);
1797 clusterInput.Fitter()->SetPrintLevel(-1+fDebugLevel);
1799 clusterInput.Fitter()->mnexcm("SET NOW", arglist, 0, ierflag);
1800 // Set starting values
1801 static Double_t vstart[6];
1802 vstart[0]=fXInit[0];
1803 vstart[1]=fYInit[0];
1804 vstart[2]=fXInit[1];
1805 vstart[3]=fYInit[1];
1806 vstart[4]=fQrInit[0];
1807 vstart[5]=fQrInit[1];
1808 // lower and upper limits
1809 static Float_t lower[6], upper[6];
1813 fSeg[1]->GetPadI(fXInit[0], fYInit[0], fZPlane, ix, iy);
1814 isec=fSeg[1]->Sector(ix, iy);
1815 dpx=fSeg[1]->Dpx(isec);
1817 fSeg[0]->GetPadI(fXInit[0], fYInit[0], fZPlane, ix, iy);
1818 isec=fSeg[0]->Sector(ix, iy);
1819 dpy=fSeg[0]->Dpy(isec);
1823 Float_t xdum, ydum, zdum;
1825 printf("\n Cluster Finder: %f %f %f %f ", fXInit[0], fXInit[1],fYInit[0], fYInit[1] );
1827 // Find save upper and lower limits
1830 for (fSeg[1]->FirstPad(fXInit[0], fYInit[0], fZPlane, dpx, 0.);
1831 fSeg[1]->MorePads(); fSeg[1]->NextPad())
1833 ix=fSeg[1]->Ix(); iy=fSeg[1]->Iy();
1834 // if (fHitMap[1]->TestHit(ix, iy) == kEmpty) continue;
1835 fSeg[1]->GetPadC(ix,iy,upper[0],ydum,zdum);
1836 if (icount ==0) lower[0]=upper[0];
1839 if (lower[0]>upper[0]) {xdum=lower[0]; lower[0]=upper[0]; upper[0]=xdum;}
1840 // vstart[0] = 0.5*(lower[0]+upper[0]);
1845 for (fSeg[0]->FirstPad(fXInit[0], fYInit[0], fZPlane, 0., dpy);
1846 fSeg[0]->MorePads(); fSeg[0]->NextPad())
1848 ix=fSeg[0]->Ix(); iy=fSeg[0]->Iy();
1849 // if (fHitMap[0]->TestHit(ix, iy) == kEmpty) continue;
1850 fSeg[0]->GetPadC(ix,iy,xdum,upper[1],zdum);
1851 if (icount ==0) lower[1]=upper[1];
1855 if (lower[1]>upper[1]) {xdum=lower[1]; lower[1]=upper[1]; upper[1]=xdum;}
1856 // vstart[1] = 0.5*(lower[1]+upper[1]);
1859 fSeg[1]->GetPadI(fXInit[1], fYInit[1], fZPlane, ix, iy);
1860 isec=fSeg[1]->Sector(ix, iy);
1861 dpx=fSeg[1]->Dpx(isec);
1862 fSeg[0]->GetPadI(fXInit[1], fYInit[1], fZPlane, ix, iy);
1863 isec=fSeg[0]->Sector(ix, iy);
1864 dpy=fSeg[0]->Dpy(isec);
1867 // Find save upper and lower limits
1871 for (fSeg[1]->FirstPad(fXInit[1], fYInit[1], fZPlane, dpx, 0);
1872 fSeg[1]->MorePads(); fSeg[1]->NextPad())
1874 ix=fSeg[1]->Ix(); iy=fSeg[1]->Iy();
1875 // if (fHitMap[1]->TestHit(ix, iy) == kEmpty) continue;
1876 fSeg[1]->GetPadC(ix,iy,upper[2],ydum,zdum);
1877 if (icount ==0) lower[2]=upper[2];
1880 if (lower[2]>upper[2]) {xdum=lower[2]; lower[2]=upper[2]; upper[2]=xdum;}
1881 // vstart[2] = 0.5*(lower[2]+upper[2]);
1885 for (fSeg[0]->FirstPad(fXInit[1], fYInit[1], fZPlane, 0, dpy);
1886 fSeg[0]-> MorePads(); fSeg[0]->NextPad())
1888 ix=fSeg[0]->Ix(); iy=fSeg[0]->Iy();
1889 // if (fHitMap[0]->TestHit(ix, iy) != kEmpty) continue;
1891 fSeg[0]->GetPadC(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;}
1898 // vstart[3] = 0.5*(lower[3]+upper[3]);
1906 static Double_t step[6]={0.0005, 0.0005, 0.0005, 0.0005, 0.001, 0.001};
1907 clusterInput.Fitter()->mnparm(0,"x1",vstart[0],step[0],lower[0],upper[0],ierflag);
1908 clusterInput.Fitter()->mnparm(1,"y1",vstart[1],step[1],lower[1],upper[1],ierflag);
1909 clusterInput.Fitter()->mnparm(2,"x2",vstart[2],step[2],lower[2],upper[2],ierflag);
1910 clusterInput.Fitter()->mnparm(3,"y2",vstart[3],step[3],lower[3],upper[3],ierflag);
1911 clusterInput.Fitter()->mnparm(4,"a0",vstart[4],step[4],lower[4],upper[4],ierflag);
1912 clusterInput.Fitter()->mnparm(5,"a1",vstart[5],step[5],lower[5],upper[5],ierflag);
1913 // ready for minimisation
1917 clusterInput.Fitter()->mnexcm("SET NOGR", arglist, 0, ierflag);
1918 clusterInput.Fitter()->mnexcm("MIGRAD", arglist, 0, ierflag);
1919 // clusterInput.Fitter()->mnexcm("EXIT" , arglist, 0, ierflag);
1920 // Get fitted parameters
1922 Double_t epxz, b1, b2;
1924 clusterInput.Fitter()->mnpout(0, chname, fXFit[0], epxz, b1, b2, ierflg);
1925 clusterInput.Fitter()->mnpout(1, chname, fYFit[0], epxz, b1, b2, ierflg);
1926 clusterInput.Fitter()->mnpout(2, chname, fXFit[1], epxz, b1, b2, ierflg);
1927 clusterInput.Fitter()->mnpout(3, chname, fYFit[1], epxz, b1, b2, ierflg);
1928 clusterInput.Fitter()->mnpout(4, chname, fQrFit[0], epxz, b1, b2, ierflg);
1929 clusterInput.Fitter()->mnpout(5, chname, fQrFit[1], epxz, b1, b2, ierflg);
1931 Double_t fmin, fedm, errdef;
1932 Int_t npari, nparx, istat;
1934 clusterInput.Fitter()->mnstat(fmin, fedm, errdef, npari, nparx, istat);
1942 void AliMUONClusterFinderVS::Split(AliMUONRawCluster* c)
1945 // One cluster for each maximum
1948 AliMUONClusterInput& clusterInput = *(AliMUONClusterInput::Instance());
1949 for (j=0; j<2; j++) {
1950 AliMUONRawCluster cnew;
1951 cnew.fGhost=c->fGhost;
1952 for (cath=0; cath<2; cath++) {
1953 cnew.fChi2[cath]=fChi2[0];
1954 // ?? why not cnew.fChi2[cath]=fChi2[cath];
1957 cnew.fNcluster[0]=-1;
1958 cnew.fNcluster[1]=fNRawClusters;
1960 cnew.fNcluster[0]=fNPeaks;
1961 cnew.fNcluster[1]=0;
1963 cnew.fMultiplicity[cath]=0;
1964 cnew.fX[cath]=Float_t(fXFit[j]);
1965 cnew.fY[cath]=Float_t(fYFit[j]);
1966 cnew.fZ[cath]=fZPlane;
1968 cnew.fQ[cath]=Int_t(clusterInput.TotalCharge(cath)*fQrFit[cath]);
1970 cnew.fQ[cath]=Int_t(clusterInput.TotalCharge(cath)*(1-fQrFit[cath]));
1972 fSeg[cath]->SetHit(fXFit[j],fYFit[j],fZPlane);
1973 for (i=0; i<fMul[cath]; i++) {
1974 cnew.fIndexMap[cnew.fMultiplicity[cath]][cath]=
1975 c->fIndexMap[i][cath];
1976 fSeg[cath]->SetPad(fIx[i][cath], fIy[i][cath]);
1977 Float_t q1=fInput->Response()->IntXY(fSeg[cath]);
1978 cnew.fContMap[i][cath]
1979 =(q1*Float_t(cnew.fQ[cath]))/Float_t(fQ[i][cath]);
1980 cnew.fMultiplicity[cath]++;
1982 FillCluster(&cnew,0,cath);
1985 cnew.fClusterType=cnew.PhysicsContribution();
1986 if (cnew.fQ[0]>0 && cnew.fQ[1]>0) AddRawCluster(cnew);
1993 // Minimisation functions
1995 void fcnS1(Int_t & /*npar*/, Double_t * /*gin*/, Double_t &f, Double_t *par, Int_t /*iflag*/)
1997 AliMUONClusterInput& clusterInput = *(AliMUONClusterInput::Instance());
2004 for (i=0; i<clusterInput.Nmul(0); i++) {
2005 Float_t q0=clusterInput.Charge(i,0);
2006 Float_t q1=clusterInput.DiscrChargeS1(i,par);
2015 void fcnCombiS1(Int_t & /*npar*/, Double_t * /*gin*/, Double_t &f, Double_t *par, Int_t /*iflag*/)
2017 AliMUONClusterInput& clusterInput = *(AliMUONClusterInput::Instance());
2024 for (cath=0; cath<2; cath++) {
2025 for (i=0; i<clusterInput.Nmul(cath); i++) {
2026 Float_t q0=clusterInput.Charge(i,cath);
2027 Float_t q1=clusterInput.DiscrChargeCombiS1(i,par,cath);
2038 void fcnS2(Int_t & /*npar*/, Double_t * /*gin*/, Double_t &f, Double_t *par, Int_t /*iflag*/)
2040 AliMUONClusterInput& clusterInput = *(AliMUONClusterInput::Instance());
2047 for (i=0; i<clusterInput.Nmul(0); i++) {
2049 Float_t q0=clusterInput.Charge(i,0);
2050 Float_t q1=clusterInput.DiscrChargeS2(i,par);
2060 void fcnCombiS2(Int_t & /*npar*/, Double_t * /*gin*/, Double_t &f, Double_t *par, Int_t /*iflag*/)
2062 AliMUONClusterInput& clusterInput = *(AliMUONClusterInput::Instance());
2068 for (cath=0; cath<2; cath++) {
2069 for (i=0; i<clusterInput.Nmul(cath); i++) {
2070 Float_t q0=clusterInput.Charge(i,cath);
2071 Float_t q1=clusterInput.DiscrChargeCombiS2(i,par,cath);
2081 void AliMUONClusterFinderVS::AddRawCluster(const AliMUONRawCluster& c)
2084 // Add a raw cluster copy to the list
2087 // AliMUON *pMUON=(AliMUON*)gAlice->GetModule("MUON");
2088 // pMUON->GetMUONData()->AddRawCluster(fInput->Chamber(),c);
2092 TClonesArray &lrawcl = *fRawClusters;
2093 new(lrawcl[fNRawClusters++]) AliMUONRawCluster(c);
2095 fprintf(stderr,"\nfNRawClusters %d\n",fNRawClusters);
2098 Bool_t AliMUONClusterFinderVS::TestTrack(Int_t t) {
2099 // Test if track was user selected
2100 if (fTrack[0]==-1 || fTrack[1]==-1) {
2102 } else if (t==fTrack[0] || t==fTrack[1]) {
2109 AliMUONClusterFinderVS& AliMUONClusterFinderVS
2110 ::operator = (const AliMUONClusterFinderVS& /*rhs*/)
2112 // Dummy assignment operator