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 **************************************************************************/
17 Revision 1.21 2001/05/18 08:54:59 morsch
18 Bug in decision on splitting corrected.
20 Revision 1.20 2001/04/12 12:22:26 morsch
21 - some numerical problems caused by pad staggering cured.
22 - treatment of 1-2 and 2-1 ghosts
23 - debuglevel > 1 prints introduced
25 Revision 1.19 2001/03/20 13:32:10 egangler
26 Code introduced to remove ghosts with the charge correlation between the 2
27 cathods. A chi2 is performed for the 2 possibilities.
28 If one gets good chi2 (with respect to the fGhostChi2Cut parameter) and the
29 other wrong chi2, the ambiguity is solved
30 If both gets good or both bad chi2, no choice is made
31 By default the fGhostChi2Cut parameter is set to solve about 70% of ghost
32 problems with about 2% errors, with the current version of the code.
35 fGhostChi2Cut is in AliMUONClusterFinderVS, with setters and getters.
36 a fDebugLevel was also introduced to switch off some of the output.
37 When an ambiguity is detected and not solved, the fGhost word in
38 AliMUONRawCluster is set to 1 or 2, depending whether both charge chi2 are
40 a DumpIndex method was also added in AliMUONRawCluster to produce a printout
44 By default, the code makes ghost check. If you want previous behaviour,
45 put in MUONrawclusters the value of SetGhostChi2Cut to infinity (1e.6) is
48 Revision 1.18 2001/01/26 21:37:53 morsch
49 Use access functions to AliMUONDigit member data.
51 Revision 1.17 2001/01/23 18:58:19 hristov
52 Initialisation of some pointers
54 Revision 1.16 2000/12/21 23:27:30 morsch
55 Error in argument list of AddRawCluster corrected.
57 Revision 1.15 2000/12/21 22:14:38 morsch
58 Clean-up of coding rule violations.
60 Revision 1.14 2000/10/23 16:03:45 morsch
61 Correct z-position of all clusters created "on the flight".
63 Revision 1.13 2000/10/23 13:38:23 morsch
64 Set correct z-coordinate when cluster is split.
66 Revision 1.12 2000/10/18 11:42:06 morsch
67 - AliMUONRawCluster contains z-position.
68 - Some clean-up of useless print statements during initialisations.
70 Revision 1.11 2000/10/06 09:04:05 morsch
71 - Dummy z-arguments in GetPadI, SetHit, FirstPad replaced by real z-coordinate
72 to make code work with slat chambers (AM)
73 - Replace GetPadI calls with unchecked x,y coordinates by pad iterator calls wherever possible.
75 Revision 1.10 2000/10/03 13:51:57 egangler
76 Removal of useless dependencies via forward declarations
78 Revision 1.9 2000/10/02 16:58:29 egangler
79 Cleaning of the code :
82 -> some useless includes removed or replaced by "class" statement
84 Revision 1.8 2000/07/03 11:54:57 morsch
85 AliMUONSegmentation and AliMUONHitMap have been replaced by AliSegmentation and AliHitMap in STEER
86 The methods GetPadIxy and GetPadXxy of AliMUONSegmentation have changed name to GetPadI and GetPadC.
88 Revision 1.7 2000/06/28 15:16:35 morsch
89 (1) Client code adapted to new method signatures in AliMUONSegmentation (see comments there)
90 to allow development of slat-muon chamber simulation and reconstruction code in the MUON
91 framework. The changes should have no side effects (mostly dummy arguments).
92 (2) Hit disintegration uses 3-dim hit coordinates to allow simulation
93 of chambers with overlapping modules (MakePadHits, Disintegration).
95 Revision 1.6 2000/06/28 12:19:18 morsch
96 More consequent seperation of global input data services (AliMUONClusterInput singleton) and the
97 cluster and hit reconstruction algorithms in AliMUONClusterFinderVS.
98 AliMUONClusterFinderVS becomes the base class for clustering and hit reconstruction.
99 It requires two cathode planes. Small modifications in the code will make it usable for
100 one cathode plane and, hence, more general (for test beam data).
101 AliMUONClusterFinder is now obsolete.
103 Revision 1.5 2000/06/28 08:06:10 morsch
104 Avoid global variables in AliMUONClusterFinderVS by seperating the input data for the fit from the
105 algorithmic part of the class. Input data resides inside the AliMUONClusterInput singleton.
106 It also naturally takes care of the TMinuit instance.
108 Revision 1.4 2000/06/27 16:18:47 gosset
109 Finally correct implementation of xm, ym, ixm, iym sizes
110 when at least three local maxima on cathode 1 or on cathode 2
112 Revision 1.3 2000/06/22 14:02:45 morsch
113 Parameterised size of xm[], ym[], ixm[], iym[] correctly implemented (PH)
114 Some HP scope problems corrected (PH)
116 Revision 1.2 2000/06/15 07:58:48 morsch
117 Code from MUON-dev joined
119 Revision 1.1.2.3 2000/06/09 21:58:33 morsch
120 Most coding rule violations corrected.
122 Revision 1.1.2.2 2000/02/15 08:33:52 morsch
123 Error in calculation of contribution map for double clusters (Split method) corrected (A.M.)
124 Error in determination of track list for double cluster (FillCluster method) corrected (A.M.)
125 Revised and extended SplitByLocalMaxima method (Isabelle Chevrot):
126 - For clusters with more than 2 maxima on one of the cathode planes all valid
127 combinations of maxima on the two cathodes are preserved. The position of the maxima is
128 taken as the hit position.
129 - New FillCluster method with 2 arguments to find tracks associated to the clusters
130 defined above added. (Method destinction by argument list not very elegant in this case,
131 should be revides (A.M.)
132 - Bug in if-statement to handle maximum 1 maximum per plane corrected
133 - Two cluster per cathode but only 1 combination valid is handled.
134 - More rigerous treatment of 1-2 and 2-1 combinations of maxima.
138 #include "AliMUONClusterFinderVS.h"
139 #include "AliMUONDigit.h"
140 #include "AliMUONRawCluster.h"
141 #include "AliSegmentation.h"
142 #include "AliMUONResponse.h"
143 #include "AliMUONClusterInput.h"
144 #include "AliMUONHitMapA1.h"
153 #include <TPostScript.h>
158 #include <Riostream.h>
160 //_____________________________________________________________________
161 // This function is minimized in the double-Mathieson fit
162 void fcnS2(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag);
163 void fcnS1(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag);
164 void fcnCombiS1(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag);
165 void fcnCombiS2(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag);
167 ClassImp(AliMUONClusterFinderVS)
169 AliMUONClusterFinderVS::AliMUONClusterFinderVS()
171 // Default constructor
172 fInput=AliMUONClusterInput::Instance();
175 fTrack[0]=fTrack[1]=-1;
176 fDebugLevel = 0; // make silent default
177 fGhostChi2Cut = 1e6; // nothing done by default
180 for(Int_t i=0; i<100; i++) {
181 for (Int_t j=0; j<2; j++) {
187 AliMUONClusterFinderVS::AliMUONClusterFinderVS(
188 const AliMUONClusterFinderVS & clusterFinder)
190 // Dummy copy Constructor
194 void AliMUONClusterFinderVS::Decluster(AliMUONRawCluster *cluster)
196 // Decluster by local maxima
197 SplitByLocalMaxima(cluster);
200 void AliMUONClusterFinderVS::SplitByLocalMaxima(AliMUONRawCluster *c)
202 // Split complex cluster by local maxima
205 fInput->SetCluster(c);
207 fMul[0]=c->fMultiplicity[0];
208 fMul[1]=c->fMultiplicity[1];
211 // dump digit information into arrays
216 for (cath=0; cath<2; cath++) {
218 for (i=0; i<fMul[cath]; i++)
221 fDig[i][cath]=fInput->Digit(cath, c->fIndexMap[i][cath]);
223 fIx[i][cath]= fDig[i][cath]->PadX();
224 fIy[i][cath]= fDig[i][cath]->PadY();
226 fQ[i][cath] = fDig[i][cath]->Signal();
227 // pad centre coordinates
229 GetPadC(fIx[i][cath], fIy[i][cath], fX[i][cath], fY[i][cath], fZ[i][cath]);
230 } // loop over cluster digits
231 } // loop over cathodes
237 // Initialise and perform mathieson fits
238 Float_t chi2, oldchi2;
239 // ++++++++++++++++++*************+++++++++++++++++++++
240 // (1) No more than one local maximum per cathode plane
241 // +++++++++++++++++++++++++++++++*************++++++++
242 if ((fNLocal[0]==1 && (fNLocal[1]==0 || fNLocal[1]==1)) ||
243 (fNLocal[0]==0 && fNLocal[1]==1)) {
244 // Perform combined single Mathieson fit
245 // Initial values for coordinates (x,y)
247 // One local maximum on cathodes 1 and 2 (X->cathode 2, Y->cathode 1)
248 if (fNLocal[0]==1 && fNLocal[1]==1) {
251 // One local maximum on cathode 1 (X,Y->cathode 1)
252 } else if (fNLocal[0]==1) {
255 // One local maximum on cathode 2 (X,Y->cathode 2)
261 fprintf(stderr,"\n cas (1) CombiSingleMathiesonFit(c)\n");
262 chi2=CombiSingleMathiesonFit(c);
263 // Int_t ndf = fgNbins[0]+fgNbins[1]-2;
264 // Float_t prob = TMath::Prob(Double_t(chi2),ndf);
265 // prob1->Fill(prob);
266 // chi2_1->Fill(chi2);
269 fprintf(stderr," chi2 %f ",chi2);
279 c->fX[0]=fSeg[0]->GetAnod(c->fX[0]);
280 c->fX[1]=fSeg[1]->GetAnod(c->fX[1]);
282 // If reasonable chi^2 add result to the list of rawclusters
285 // If not try combined double Mathieson Fit
287 fprintf(stderr," MAUVAIS CHI2 !!!\n");
288 if (fNLocal[0]==1 && fNLocal[1]==1) {
289 fXInit[0]=fX[fIndLocal[0][1]][1];
290 fYInit[0]=fY[fIndLocal[0][0]][0];
291 fXInit[1]=fX[fIndLocal[0][1]][1];
292 fYInit[1]=fY[fIndLocal[0][0]][0];
293 } else if (fNLocal[0]==1) {
294 fXInit[0]=fX[fIndLocal[0][0]][0];
295 fYInit[0]=fY[fIndLocal[0][0]][0];
296 fXInit[1]=fX[fIndLocal[0][0]][0];
297 fYInit[1]=fY[fIndLocal[0][0]][0];
299 fXInit[0]=fX[fIndLocal[0][1]][1];
300 fYInit[0]=fY[fIndLocal[0][1]][1];
301 fXInit[1]=fX[fIndLocal[0][1]][1];
302 fYInit[1]=fY[fIndLocal[0][1]][1];
305 // Initial value for charge ratios
309 fprintf(stderr,"\n cas (1) CombiDoubleMathiesonFit(c)\n");
310 chi2=CombiDoubleMathiesonFit(c);
311 // Int_t ndf = fgNbins[0]+fgNbins[1]-6;
312 // Float_t prob = TMath::Prob(chi2,ndf);
313 // prob2->Fill(prob);
314 // chi2_2->Fill(chi2);
316 // Was this any better ??
317 fprintf(stderr," Old and new chi2 %f %f ", oldchi2, chi2);
318 if (fFitStat!=0 && chi2>0 && (2.*chi2 < oldchi2)) {
319 fprintf(stderr," Split\n");
320 // Split cluster into two according to fit result
323 fprintf(stderr," Don't Split\n");
329 // +++++++++++++++++++++++++++++++++++++++
330 // (2) Two local maxima per cathode plane
331 // +++++++++++++++++++++++++++++++++++++++
332 } else if (fNLocal[0]==2 && fNLocal[1]==2) {
334 // Let's look for ghosts first
336 Float_t xm[4][2], ym[4][2];
337 Float_t dpx, dpy, dx, dy;
338 Int_t ixm[4][2], iym[4][2];
339 Int_t isec, im1, im2, ico;
341 // Form the 2x2 combinations
342 // 0-0, 0-1, 1-0, 1-1
344 for (im1=0; im1<2; im1++) {
345 for (im2=0; im2<2; im2++) {
346 xm[ico][0]=fX[fIndLocal[im1][0]][0];
347 ym[ico][0]=fY[fIndLocal[im1][0]][0];
348 xm[ico][1]=fX[fIndLocal[im2][1]][1];
349 ym[ico][1]=fY[fIndLocal[im2][1]][1];
351 ixm[ico][0]=fIx[fIndLocal[im1][0]][0];
352 iym[ico][0]=fIy[fIndLocal[im1][0]][0];
353 ixm[ico][1]=fIx[fIndLocal[im2][1]][1];
354 iym[ico][1]=fIy[fIndLocal[im2][1]][1];
358 // ico = 0 : first local maximum on cathodes 1 and 2
359 // ico = 1 : fisrt local maximum on cathode 1 and second on cathode 2
360 // ico = 2 : second local maximum on cathode 1 and first on cathode 1
361 // ico = 3 : second local maximum on cathodes 1 and 2
363 // Analyse the combinations and keep those that are possible !
364 // For each combination check consistency in x and y
367 Float_t dr[4] = {1.e4, 1.e4, 1.e4, 1.e4};
370 // In case of staggering maxima are displaced by exactly half the pad-size in y.
371 // We have to take into account the numerical precision in the consistency check;
374 for (ico=0; ico<4; ico++) {
375 accepted[ico]=kFALSE;
376 // cathode one: x-coordinate
377 isec=fSeg[0]->Sector(ixm[ico][0], iym[ico][0]);
378 dpx=fSeg[0]->Dpx(isec)/2.;
379 dx=TMath::Abs(xm[ico][0]-xm[ico][1]);
380 // cathode two: y-coordinate
381 isec=fSeg[1]->Sector(ixm[ico][1], iym[ico][1]);
382 dpy=fSeg[1]->Dpy(isec)/2.;
383 dy=TMath::Abs(ym[ico][0]-ym[ico][1]);
385 printf("\n %i %f %f %f %f %f %f \n", ico, ym[ico][0], ym[ico][1], dy, dpy, dx, dpx );
386 if ((dx <= dpx) && (dy <= dpy+eps)) {
389 dr[ico] = TMath::Sqrt(dx*dx+dy*dy);
393 accepted[ico]=kFALSE;
396 printf("\n iacc= %d:\n", iacc);
398 if (accepted[0] && accepted[1]) {
399 if (dr[0] >= dr[1]) {
406 if (accepted[2] && accepted[3]) {
407 if (dr[2] >= dr[3]) {
414 // eliminate one candidate
418 for (ico=0; ico<4; ico++) {
419 if (accepted[ico] && dr[ico] > drmax) {
425 accepted[icobad] = kFALSE;
431 printf("\n iacc= %d:\n", iacc);
434 fprintf(stderr,"\n iacc=2: No problem ! \n");
435 } else if (iacc==4) {
436 fprintf(stderr,"\n iacc=4: Ok, but ghost problem !!! \n");
437 } else if (iacc==0) {
438 fprintf(stderr,"\n iacc=0: I don't know what to do with this !!!!!!!!! \n");
442 // Initial value for charge ratios
443 fQrInit[0]=Float_t(fQ[fIndLocal[0][0]][0])/
444 Float_t(fQ[fIndLocal[0][0]][0]+fQ[fIndLocal[1][0]][0]);
445 fQrInit[1]=Float_t(fQ[fIndLocal[0][1]][1])/
446 Float_t(fQ[fIndLocal[0][1]][1]+fQ[fIndLocal[1][1]][1]);
448 // ******* iacc = 0 *******
449 // No combinations found between the 2 cathodes
450 // We keep the center of gravity of the cluster
455 // ******* iacc = 1 *******
456 // Only one combination found between the 2 cathodes
458 // Initial values for the 2 maxima (x,y)
460 // 1 maximum is initialised with the maximum of the combination found (X->cathode 2, Y->cathode 1)
461 // 1 maximum is initialised with the other maximum of the first cathode
463 fprintf(stderr,"ico=0\n");
468 } else if (accepted[1]){
469 fprintf(stderr,"ico=1\n");
474 } else if (accepted[2]){
475 fprintf(stderr,"ico=2\n");
480 } else if (accepted[3]){
481 fprintf(stderr,"ico=3\n");
488 fprintf(stderr,"\n cas (2) CombiDoubleMathiesonFit(c)\n");
489 chi2=CombiDoubleMathiesonFit(c);
490 // Int_t ndf = fgNbins[0]+fgNbins[1]-6;
491 // Float_t prob = TMath::Prob(chi2,ndf);
492 // prob2->Fill(prob);
493 // chi2_2->Fill(chi2);
495 fprintf(stderr," chi2 %f\n",chi2);
497 // If reasonable chi^2 add result to the list of rawclusters
502 // 1 maximum is initialised with the maximum of the combination found (X->cathode 2, Y->cathode 1)
503 // 1 maximum is initialised with the other maximum of the second cathode
505 fprintf(stderr,"ico=0\n");
510 } else if (accepted[1]){
511 fprintf(stderr,"ico=1\n");
516 } else if (accepted[2]){
517 fprintf(stderr,"ico=2\n");
522 } else if (accepted[3]){
523 fprintf(stderr,"ico=3\n");
530 fprintf(stderr,"\n cas (2) CombiDoubleMathiesonFit(c)\n");
531 chi2=CombiDoubleMathiesonFit(c);
532 // Int_t ndf = fgNbins[0]+fgNbins[1]-6;
533 // Float_t prob = TMath::Prob(chi2,ndf);
534 // prob2->Fill(prob);
535 // chi2_2->Fill(chi2);
537 fprintf(stderr," chi2 %f\n",chi2);
539 // If reasonable chi^2 add result to the list of rawclusters
543 //We keep only the combination found (X->cathode 2, Y->cathode 1)
544 for (Int_t ico=0; ico<2; ico++) {
546 AliMUONRawCluster cnew;
548 for (cath=0; cath<2; cath++) {
549 cnew.fX[cath]=Float_t(xm[ico][1]);
550 cnew.fY[cath]=Float_t(ym[ico][0]);
551 cnew.fZ[cath]=fZPlane;
553 cnew.fMultiplicity[cath]=c->fMultiplicity[cath];
554 for (i=0; i<fMul[cath]; i++) {
555 cnew.fIndexMap[i][cath]=c->fIndexMap[i][cath];
556 fSeg[cath]->SetPad(fIx[i][cath], fIy[i][cath]);
558 fprintf(stderr,"\nRawCluster %d cath %d\n",ico,cath);
559 fprintf(stderr,"mult_av %d\n",c->fMultiplicity[cath]);
560 FillCluster(&cnew,cath);
562 cnew.fClusterType=cnew.PhysicsContribution();
571 // ******* iacc = 2 *******
572 // Two combinations found between the 2 cathodes
574 // Was the same maximum taken twice
575 if ((accepted[0]&&accepted[1]) || (accepted[2]&&accepted[3])) {
576 fprintf(stderr,"\n Maximum taken twice !!!\n");
578 // Have a try !! with that
579 if (accepted[0]&&accepted[3]) {
591 fprintf(stderr,"\n cas (2) CombiDoubleMathiesonFit(c)\n");
592 chi2=CombiDoubleMathiesonFit(c);
593 // Int_t ndf = fgNbins[0]+fgNbins[1]-6;
594 // Float_t prob = TMath::Prob(chi2,ndf);
595 // prob2->Fill(prob);
596 // chi2_2->Fill(chi2);
600 // No ghosts ! No Problems ! - Perform one fit only !
601 if (accepted[0]&&accepted[3]) {
613 fprintf(stderr,"\n cas (2) CombiDoubleMathiesonFit(c)\n");
614 chi2=CombiDoubleMathiesonFit(c);
615 // Int_t ndf = fgNbins[0]+fgNbins[1]-6;
616 // Float_t prob = TMath::Prob(chi2,ndf);
617 // prob2->Fill(prob);
618 // chi2_2->Fill(chi2);
620 fprintf(stderr," chi2 %f\n",chi2);
624 // ******* iacc = 4 *******
625 // Four combinations found between the 2 cathodes
627 } else if (iacc==4) {
628 // Perform fits for the two possibilities !!
629 // Accept if charges are compatible on both cathodes
630 // If none are compatible, keep everything
636 fprintf(stderr,"\n cas (2) CombiDoubleMathiesonFit(c)\n");
637 chi2=CombiDoubleMathiesonFit(c);
638 // Int_t ndf = fgNbins[0]+fgNbins[1]-6;
639 // Float_t prob = TMath::Prob(chi2,ndf);
640 // prob2->Fill(prob);
641 // chi2_2->Fill(chi2);
643 fprintf(stderr," chi2 %f\n",chi2);
644 // store results of fit and postpone decision
645 Double_t sXFit[2],sYFit[2],sQrFit[2];
647 for (Int_t i=0;i<2;i++) {
658 fprintf(stderr,"\n cas (2) CombiDoubleMathiesonFit(c)\n");
659 chi2=CombiDoubleMathiesonFit(c);
660 // ndf = fgNbins[0]+fgNbins[1]-6;
661 // prob = TMath::Prob(chi2,ndf);
662 // prob2->Fill(prob);
663 // chi2_2->Fill(chi2);
665 fprintf(stderr," chi2 %f\n",chi2);
666 // We have all informations to perform the decision
667 // Compute the chi2 for the 2 possibilities
668 Float_t chi2fi,chi2si,chi2f,chi2s;
670 chi2f = (TMath::Log(fInput->TotalCharge(0)*fQrFit[0]
671 / (fInput->TotalCharge(1)*fQrFit[1]) )
672 / fInput->Response()->ChargeCorrel() );
674 chi2fi = (TMath::Log(fInput->TotalCharge(0)*(1-fQrFit[0])
675 / (fInput->TotalCharge(1)*(1-fQrFit[1])) )
676 / fInput->Response()->ChargeCorrel() );
677 chi2f += chi2fi*chi2fi;
679 chi2s = (TMath::Log(fInput->TotalCharge(0)*sQrFit[0]
680 / (fInput->TotalCharge(1)*sQrFit[1]) )
681 / fInput->Response()->ChargeCorrel() );
683 chi2si = (TMath::Log(fInput->TotalCharge(0)*(1-sQrFit[0])
684 / (fInput->TotalCharge(1)*(1-sQrFit[1])) )
685 / fInput->Response()->ChargeCorrel() );
686 chi2s += chi2si*chi2si;
688 // usefull to store the charge matching chi2 in the cluster
689 // fChi2[0]=sChi2[1]=chi2f;
690 // fChi2[1]=sChi2[0]=chi2s;
692 if (chi2f<=fGhostChi2Cut && chi2s<=fGhostChi2Cut)
694 if (chi2f>fGhostChi2Cut && chi2s>fGhostChi2Cut) {
700 if (chi2f<=fGhostChi2Cut)
702 if (chi2s<=fGhostChi2Cut) {
703 // retreive saved values
704 for (Int_t i=0;i<2;i++) {
715 } else if (fNLocal[0]==2 && fNLocal[1]==1) {
716 // ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
717 // (3) Two local maxima on cathode 1 and one maximum on cathode 2
718 // ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
720 Float_t xm[4][2], ym[4][2];
721 Float_t dpx, dpy, dx, dy;
722 Int_t ixm[4][2], iym[4][2];
723 Int_t isec, im1, ico;
725 // Form the 2x2 combinations
726 // 0-0, 0-1, 1-0, 1-1
728 for (im1=0; im1<2; im1++) {
729 xm[ico][0]=fX[fIndLocal[im1][0]][0];
730 ym[ico][0]=fY[fIndLocal[im1][0]][0];
731 xm[ico][1]=fX[fIndLocal[0][1]][1];
732 ym[ico][1]=fY[fIndLocal[0][1]][1];
734 ixm[ico][0]=fIx[fIndLocal[im1][0]][0];
735 iym[ico][0]=fIy[fIndLocal[im1][0]][0];
736 ixm[ico][1]=fIx[fIndLocal[0][1]][1];
737 iym[ico][1]=fIy[fIndLocal[0][1]][1];
740 // ico = 0 : first local maximum on cathodes 1 and 2
741 // ico = 1 : second local maximum on cathode 1 and first on cathode 2
743 // Analyse the combinations and keep those that are possible !
744 // For each combination check consistency in x and y
748 // In case of staggering maxima are displaced by exactly half the pad-size in y.
749 // We have to take into account the numerical precision in the consistency check;
753 for (ico=0; ico<2; ico++) {
754 accepted[ico]=kFALSE;
755 isec=fSeg[0]->Sector(ixm[ico][0], iym[ico][0]);
756 dpx=fSeg[0]->Dpx(isec)/2.;
757 dx=TMath::Abs(xm[ico][0]-xm[ico][1]);
758 isec=fSeg[1]->Sector(ixm[ico][1], iym[ico][1]);
759 dpy=fSeg[1]->Dpy(isec)/2.;
760 dy=TMath::Abs(ym[ico][0]-ym[ico][1]);
762 printf("\n %i %f %f %f %f \n", ico, ym[ico][0], ym[ico][1], dy, dpy );
763 if ((dx <= dpx) && (dy <= dpy+eps)) {
769 accepted[ico]=kFALSE;
777 // Initial value for charge ratios
778 fQrInit[0]=Float_t(fQ[fIndLocal[0][0]][0])/
779 Float_t(fQ[fIndLocal[0][0]][0]+fQ[fIndLocal[1][0]][0]);
780 fQrInit[1]=fQrInit[0];
782 if (accepted[0] && accepted[1]) {
784 fXInit[0]=0.5*(xm[0][1]+xm[0][0]);
786 fXInit[1]=0.5*(xm[0][1]+xm[1][0]);
790 chi23=CombiDoubleMathiesonFit(c);
799 } else if (accepted[0]) {
804 chi21=CombiDoubleMathiesonFit(c);
805 // Int_t ndf = fgNbins[0]+fgNbins[1]-6;
806 // Float_t prob = TMath::Prob(chi2,ndf);
807 // prob2->Fill(prob);
808 // chi2_2->Fill(chi21);
810 fprintf(stderr," chi2 %f\n",chi21);
811 if (chi21<10) Split(c);
812 } else if (accepted[1]) {
817 chi22=CombiDoubleMathiesonFit(c);
818 // Int_t ndf = fgNbins[0]+fgNbins[1]-6;
819 // Float_t prob = TMath::Prob(chi2,ndf);
820 // prob2->Fill(prob);
821 // chi2_2->Fill(chi22);
823 fprintf(stderr," chi2 %f\n",chi22);
824 if (chi22<10) Split(c);
827 if (chi21 > 10 && chi22 > 10 && chi23 > 10) {
828 // We keep only the combination found (X->cathode 2, Y->cathode 1)
829 for (Int_t ico=0; ico<2; ico++) {
831 AliMUONRawCluster cnew;
833 for (cath=0; cath<2; cath++) {
834 cnew.fX[cath]=Float_t(xm[ico][1]);
835 cnew.fY[cath]=Float_t(ym[ico][0]);
836 cnew.fZ[cath]=fZPlane;
837 cnew.fMultiplicity[cath]=c->fMultiplicity[cath];
838 for (i=0; i<fMul[cath]; i++) {
839 cnew.fIndexMap[i][cath]=c->fIndexMap[i][cath];
840 fSeg[cath]->SetPad(fIx[i][cath], fIy[i][cath]);
842 fprintf(stderr,"\nRawCluster %d cath %d\n",ico,cath);
843 fprintf(stderr,"mult_av %d\n",c->fMultiplicity[cath]);
844 FillCluster(&cnew,cath);
846 cnew.fClusterType=cnew.PhysicsContribution();
853 // +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
854 // (3') One local maximum on cathode 1 and two maxima on cathode 2
855 // +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
856 } else if (fNLocal[0]==1 && fNLocal[1]==2) {
857 Float_t xm[4][2], ym[4][2];
858 Float_t dpx, dpy, dx, dy;
859 Int_t ixm[4][2], iym[4][2];
860 Int_t isec, im1, ico;
862 // Form the 2x2 combinations
863 // 0-0, 0-1, 1-0, 1-1
865 for (im1=0; im1<2; im1++) {
866 xm[ico][0]=fX[fIndLocal[0][0]][0];
867 ym[ico][0]=fY[fIndLocal[0][0]][0];
868 xm[ico][1]=fX[fIndLocal[im1][1]][1];
869 ym[ico][1]=fY[fIndLocal[im1][1]][1];
871 ixm[ico][0]=fIx[fIndLocal[0][0]][0];
872 iym[ico][0]=fIy[fIndLocal[0][0]][0];
873 ixm[ico][1]=fIx[fIndLocal[im1][1]][1];
874 iym[ico][1]=fIy[fIndLocal[im1][1]][1];
877 // ico = 0 : first local maximum on cathodes 1 and 2
878 // ico = 1 : first local maximum on cathode 1 and second on cathode 2
880 // Analyse the combinations and keep those that are possible !
881 // For each combination check consistency in x and y
885 // In case of staggering maxima are displaced by exactly half the pad-size in y.
886 // We have to take into account the numerical precision in the consistency check;
890 for (ico=0; ico<2; ico++) {
891 accepted[ico]=kFALSE;
892 isec=fSeg[0]->Sector(ixm[ico][0], iym[ico][0]);
893 dpx=fSeg[0]->Dpx(isec)/2.;
894 dx=TMath::Abs(xm[ico][0]-xm[ico][1]);
895 isec=fSeg[1]->Sector(ixm[ico][1], iym[ico][1]);
896 dpy=fSeg[1]->Dpy(isec)/2.;
897 dy=TMath::Abs(ym[ico][0]-ym[ico][1]);
899 printf("\n %i %f %f %f %f \n", ico, ym[ico][0], ym[ico][1], dy, dpy );
900 if ((dx <= dpx) && (dy <= dpy+eps)) {
903 fprintf(stderr,"ico %d\n",ico);
907 accepted[ico]=kFALSE;
915 fQrInit[1]=Float_t(fQ[fIndLocal[0][1]][1])/
916 Float_t(fQ[fIndLocal[0][1]][1]+fQ[fIndLocal[1][1]][1]);
918 fQrInit[0]=fQrInit[1];
921 if (accepted[0] && accepted[1]) {
923 fYInit[0]=0.5*(ym[0][0]+ym[0][1]);
925 fYInit[1]=0.5*(ym[0][0]+ym[1][1]);
928 chi23=CombiDoubleMathiesonFit(c);
937 } else if (accepted[0]) {
942 chi21=CombiDoubleMathiesonFit(c);
943 // Int_t ndf = fgNbins[0]+fgNbins[1]-6;
944 // Float_t prob = TMath::Prob(chi2,ndf);
945 // prob2->Fill(prob);
946 // chi2_2->Fill(chi21);
948 fprintf(stderr," chi2 %f\n",chi21);
949 if (chi21<10) Split(c);
950 } else if (accepted[1]) {
955 chi22=CombiDoubleMathiesonFit(c);
956 // Int_t ndf = fgNbins[0]+fgNbins[1]-6;
957 // Float_t prob = TMath::Prob(chi2,ndf);
958 // prob2->Fill(prob);
959 // chi2_2->Fill(chi22);
961 fprintf(stderr," chi2 %f\n",chi22);
962 if (chi22<10) Split(c);
965 if (chi21 > 10 && chi22 > 10 && chi23 > 10) {
966 //We keep only the combination found (X->cathode 2, Y->cathode 1)
967 for (Int_t ico=0; ico<2; ico++) {
969 AliMUONRawCluster cnew;
971 for (cath=0; cath<2; cath++) {
972 cnew.fX[cath]=Float_t(xm[ico][1]);
973 cnew.fY[cath]=Float_t(ym[ico][0]);
974 cnew.fZ[cath]=fZPlane;
975 cnew.fMultiplicity[cath]=c->fMultiplicity[cath];
976 for (i=0; i<fMul[cath]; i++) {
977 cnew.fIndexMap[i][cath]=c->fIndexMap[i][cath];
978 fSeg[cath]->SetPad(fIx[i][cath], fIy[i][cath]);
980 fprintf(stderr,"\nRawCluster %d cath %d\n",ico,cath);
981 fprintf(stderr,"mult_av %d\n",c->fMultiplicity[cath]);
982 FillCluster(&cnew,cath);
984 cnew.fClusterType=cnew.PhysicsContribution();
991 // ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
992 // (4) At least three local maxima on cathode 1 or on cathode 2
993 // ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
994 } else if (fNLocal[0]>2 || fNLocal[1]>2) {
995 Int_t param = fNLocal[0]*fNLocal[1];
998 Float_t ** xm = new Float_t * [param];
999 for (ii=0; ii<param; ii++) xm[ii]=new Float_t [2];
1000 Float_t ** ym = new Float_t * [param];
1001 for (ii=0; ii<param; ii++) ym[ii]=new Float_t [2];
1002 Int_t ** ixm = new Int_t * [param];
1003 for (ii=0; ii<param; ii++) ixm[ii]=new Int_t [2];
1004 Int_t ** iym = new Int_t * [param];
1005 for (ii=0; ii<param; ii++) iym[ii]=new Int_t [2];
1008 Float_t dpx, dpy, dx, dy;
1011 for (Int_t im1=0; im1<fNLocal[0]; im1++) {
1012 for (Int_t im2=0; im2<fNLocal[1]; im2++) {
1013 xm[ico][0]=fX[fIndLocal[im1][0]][0];
1014 ym[ico][0]=fY[fIndLocal[im1][0]][0];
1015 xm[ico][1]=fX[fIndLocal[im2][1]][1];
1016 ym[ico][1]=fY[fIndLocal[im2][1]][1];
1018 ixm[ico][0]=fIx[fIndLocal[im1][0]][0];
1019 iym[ico][0]=fIy[fIndLocal[im1][0]][0];
1020 ixm[ico][1]=fIx[fIndLocal[im2][1]][1];
1021 iym[ico][1]=fIy[fIndLocal[im2][1]][1];
1028 fprintf(stderr,"nIco %d\n",nIco);
1029 for (ico=0; ico<nIco; ico++) {
1031 fprintf(stderr,"ico = %d\n",ico);
1032 isec=fSeg[0]->Sector(ixm[ico][0], iym[ico][0]);
1033 dpx=fSeg[0]->Dpx(isec)/2.;
1034 dx=TMath::Abs(xm[ico][0]-xm[ico][1]);
1035 isec=fSeg[1]->Sector(ixm[ico][1], iym[ico][1]);
1036 dpy=fSeg[1]->Dpy(isec)/2.;
1037 dy=TMath::Abs(ym[ico][0]-ym[ico][1]);
1039 fprintf(stderr,"dx %f dpx %f dy %f dpy %f\n",dx,dpx,dy,dpy);
1040 fprintf(stderr," X %f Y %f\n",xm[ico][1],ym[ico][0]);
1042 if ((dx <= dpx) && (dy <= dpy)) {
1044 fprintf(stderr,"ok\n");
1046 AliMUONRawCluster cnew;
1047 for (cath=0; cath<2; cath++) {
1048 cnew.fX[cath]=Float_t(xm[ico][1]);
1049 cnew.fY[cath]=Float_t(ym[ico][0]);
1050 cnew.fZ[cath]=fZPlane;
1051 cnew.fMultiplicity[cath]=c->fMultiplicity[cath];
1052 for (i=0; i<fMul[cath]; i++) {
1053 cnew.fIndexMap[i][cath]=c->fIndexMap[i][cath];
1054 fSeg[cath]->SetPad(fIx[i][cath], fIy[i][cath]);
1056 FillCluster(&cnew,cath);
1058 cnew.fClusterType=cnew.PhysicsContribution();
1059 AddRawCluster(cnew);
1070 void AliMUONClusterFinderVS::FindLocalMaxima(AliMUONRawCluster* c)
1072 // Find all local maxima of a cluster
1074 printf("\n Find Local maxima !");
1078 Int_t cath, cath1; // loops over cathodes
1079 Int_t i; // loops over digits
1080 Int_t j; // loops over cathodes
1082 // Find local maxima
1084 // counters for number of local maxima
1085 fNLocal[0]=fNLocal[1]=0;
1086 // flags digits as local maximum
1087 Bool_t isLocal[100][2];
1088 for (i=0; i<100;i++) {
1089 isLocal[i][0]=isLocal[i][1]=kFALSE;
1091 // number of next neighbours and arrays to store them
1094 // loop over cathodes
1095 for (cath=0; cath<2; cath++) {
1096 // loop over cluster digits
1097 for (i=0; i<fMul[cath]; i++) {
1098 // get neighbours for that digit and assume that it is local maximum
1099 fSeg[cath]->Neighbours(fIx[i][cath], fIy[i][cath], &nn, x, y);
1100 isLocal[i][cath]=kTRUE;
1101 Int_t isec= fSeg[cath]->Sector(fIx[i][cath], fIy[i][cath]);
1102 Float_t a0 = fSeg[cath]->Dpx(isec)*fSeg[cath]->Dpy(isec);
1103 // loop over next neighbours, if at least one neighbour has higher charger assumption
1104 // digit is not local maximum
1105 for (j=0; j<nn; j++) {
1106 if (fHitMap[cath]->TestHit(x[j], y[j])==kEmpty) continue;
1107 digt=(AliMUONDigit*) fHitMap[cath]->GetHit(x[j], y[j]);
1108 isec=fSeg[cath]->Sector(x[j], y[j]);
1109 Float_t a1 = fSeg[cath]->Dpx(isec)*fSeg[cath]->Dpy(isec);
1110 if (digt->Signal()/a1 > fQ[i][cath]/a0) {
1111 isLocal[i][cath]=kFALSE;
1114 // handle special case of neighbouring pads with equal signal
1115 } else if (digt->Signal() == fQ[i][cath]) {
1116 if (fNLocal[cath]>0) {
1117 for (Int_t k=0; k<fNLocal[cath]; k++) {
1118 if (x[j]==fIx[fIndLocal[k][cath]][cath]
1119 && y[j]==fIy[fIndLocal[k][cath]][cath])
1121 isLocal[i][cath]=kFALSE;
1123 } // loop over local maxima
1124 } // are there already local maxima
1126 } // loop over next neighbours
1127 if (isLocal[i][cath]) {
1128 fIndLocal[fNLocal[cath]][cath]=i;
1131 } // loop over all digits
1132 } // loop over cathodes
1135 printf("\n Found %d %d %d %d local Maxima\n",
1136 fNLocal[0], fNLocal[1], fMul[0], fMul[1]);
1137 fprintf(stderr,"\n Cathode 1 local Maxima %d Multiplicite %d\n",fNLocal[0], fMul[0]);
1138 fprintf(stderr," Cathode 2 local Maxima %d Multiplicite %d\n",fNLocal[1], fMul[1]);
1144 if (fNLocal[1]==2 && (fNLocal[0]==1 || fNLocal[0]==0)) {
1145 Int_t iback=fNLocal[0];
1147 // Two local maxima on cathode 2 and one maximum on cathode 1
1148 // Look for local maxima considering up and down neighbours on the 1st cathode only
1150 // Loop over cluster digits
1154 for (i=0; i<fMul[cath]; i++) {
1155 isec=fSeg[cath]->Sector(fIx[i][cath],fIy[i][cath]);
1156 dpy=fSeg[cath]->Dpy(isec);
1157 dpx=fSeg[cath]->Dpx(isec);
1158 if (isLocal[i][cath]) continue;
1159 // Pad position should be consistent with position of local maxima on the opposite cathode
1160 if ((TMath::Abs(fX[i][cath]-fX[fIndLocal[0][cath1]][cath1]) > dpx/2.) &&
1161 (TMath::Abs(fX[i][cath]-fX[fIndLocal[1][cath1]][cath1]) > dpx/2.))
1164 // get neighbours for that digit and assume that it is local maximum
1165 isLocal[i][cath]=kTRUE;
1166 // compare signal to that on the two neighbours on the left and on the right
1167 // iNN counts the number of neighbours with signal, it should be 1 or 2
1171 ->FirstPad(fX[i][cath], fY[i][cath], fZPlane, 0., dpy);
1177 ix = fSeg[cath]->Ix();
1178 iy = fSeg[cath]->Iy();
1179 // skip the current pad
1180 if (iy == fIy[i][cath]) continue;
1182 if (fHitMap[cath]->TestHit(ix, iy)!=kEmpty) {
1184 digt=(AliMUONDigit*) fHitMap[cath]->GetHit(ix,iy);
1185 if (digt->Signal() > fQ[i][cath]) isLocal[i][cath]=kFALSE;
1187 } // Loop over pad neighbours in y
1188 if (isLocal[i][cath] && iNN>0) {
1189 fIndLocal[fNLocal[cath]][cath]=i;
1192 } // loop over all digits
1193 // if one additional maximum has been found we are happy
1194 // if more maxima have been found restore the previous situation
1197 "\n New search gives %d local maxima for cathode 1 \n",
1200 " %d local maxima for cathode 2 \n",
1203 if (fNLocal[cath]>2) {
1204 fNLocal[cath]=iback;
1207 } // 1,2 local maxima
1209 if (fNLocal[0]==2 && (fNLocal[1]==1 || fNLocal[1]==0)) {
1210 Int_t iback=fNLocal[1];
1212 // Two local maxima on cathode 1 and one maximum on cathode 2
1213 // Look for local maxima considering left and right neighbours on the 2nd cathode only
1216 Float_t eps = 1.e-5;
1219 // Loop over cluster digits
1220 for (i=0; i<fMul[cath]; i++) {
1221 isec=fSeg[cath]->Sector(fIx[i][cath],fIy[i][cath]);
1222 dpx=fSeg[cath]->Dpx(isec);
1223 dpy=fSeg[cath]->Dpy(isec);
1224 if (isLocal[i][cath]) continue;
1225 // Pad position should be consistent with position of local maxima on the opposite cathode
1226 if ((TMath::Abs(fY[i][cath]-fY[fIndLocal[0][cath1]][cath1]) > dpy/2.+eps) &&
1227 (TMath::Abs(fY[i][cath]-fY[fIndLocal[1][cath1]][cath1]) > dpy/2.+eps))
1231 // get neighbours for that digit and assume that it is local maximum
1232 isLocal[i][cath]=kTRUE;
1233 // compare signal to that on the two neighbours on the left and on the right
1235 // iNN counts the number of neighbours with signal, it should be 1 or 2
1238 ->FirstPad(fX[i][cath], fY[i][cath], fZPlane, dpx, 0.);
1245 ix = fSeg[cath]->Ix();
1246 iy = fSeg[cath]->Iy();
1248 // skip the current pad
1249 if (ix == fIx[i][cath]) continue;
1251 if (fHitMap[cath]->TestHit(ix, iy)!=kEmpty) {
1253 digt=(AliMUONDigit*) fHitMap[cath]->GetHit(ix,iy);
1254 if (digt->Signal() > fQ[i][cath]) isLocal[i][cath]=kFALSE;
1256 } // Loop over pad neighbours in x
1257 if (isLocal[i][cath] && iNN>0) {
1258 fIndLocal[fNLocal[cath]][cath]=i;
1261 } // loop over all digits
1262 // if one additional maximum has been found we are happy
1263 // if more maxima have been found restore the previous situation
1265 fprintf(stderr,"\n New search gives %d local maxima for cathode 1 \n",fNLocal[0]);
1266 fprintf(stderr,"\n %d local maxima for cathode 2 \n",fNLocal[1]);
1267 printf("\n New search gives %d %d \n",fNLocal[0],fNLocal[1]);
1269 if (fNLocal[cath]>2) {
1270 fNLocal[cath]=iback;
1272 } // 2,1 local maxima
1276 void AliMUONClusterFinderVS::FillCluster(AliMUONRawCluster* c, Int_t flag, Int_t cath)
1279 // Completes cluster information starting from list of digits
1286 c->fPeakSignal[cath]=c->fPeakSignal[0];
1288 c->fPeakSignal[cath]=0;
1299 fprintf(stderr,"\n fPeakSignal %d\n",c->fPeakSignal[cath]);
1300 for (Int_t i=0; i<c->fMultiplicity[cath]; i++)
1302 dig= fInput->Digit(cath,c->fIndexMap[i][cath]);
1303 ix=dig->PadX()+c->fOffsetMap[i][cath];
1305 Int_t q=dig->Signal();
1306 if (!flag) q=Int_t(q*c->fContMap[i][cath]);
1307 // fprintf(stderr,"q %d c->fPeakSignal[ %d ] %d\n",q,cath,c->fPeakSignal[cath]);
1308 if (dig->Physics() >= dig->Signal()) {
1309 c->fPhysicsMap[i]=2;
1310 } else if (dig->Physics() == 0) {
1311 c->fPhysicsMap[i]=0;
1312 } else c->fPhysicsMap[i]=1;
1316 fprintf(stderr,"q %d c->fPeakSignal[cath] %d\n",q,c->fPeakSignal[cath]);
1317 // peak signal and track list
1318 if (q>c->fPeakSignal[cath]) {
1319 c->fPeakSignal[cath]=q;
1320 c->fTracks[0]=dig->Hit();
1321 c->fTracks[1]=dig->Track(0);
1322 c->fTracks[2]=dig->Track(1);
1323 // fprintf(stderr," c->fTracks[0] %d c->fTracks[1] %d\n",dig->fHit,dig->fTracks[0]);
1327 fSeg[cath]->GetPadC(ix, iy, x, y, z);
1332 } // loop over digits
1334 fprintf(stderr," fin du cluster c\n");
1338 c->fX[cath]/=c->fQ[cath];
1340 c->fX[cath]=fSeg[cath]->GetAnod(c->fX[cath]);
1341 c->fY[cath]/=c->fQ[cath];
1343 // apply correction to the coordinate along the anode wire
1347 fSeg[cath]->GetPadI(x, y, fZPlane, ix, iy);
1348 fSeg[cath]->GetPadC(ix, iy, x, y, z);
1349 Int_t isec=fSeg[cath]->Sector(ix,iy);
1350 TF1* cogCorr = fSeg[cath]->CorrFunc(isec-1);
1353 Float_t yOnPad=(c->fY[cath]-y)/fSeg[cath]->Dpy(isec);
1354 c->fY[cath]=c->fY[cath]-cogCorr->Eval(yOnPad, 0, 0);
1359 void AliMUONClusterFinderVS::FillCluster(AliMUONRawCluster* c, Int_t cath)
1362 // Completes cluster information starting from list of digits
1372 Float_t xpad, ypad, zpad;
1375 for (Int_t i=0; i<c->fMultiplicity[cath]; i++)
1377 dig = fInput->Digit(cath,c->fIndexMap[i][cath]);
1379 GetPadC(dig->PadX(),dig->PadY(),xpad,ypad, zpad);
1381 fprintf(stderr,"x %f y %f cx %f cy %f\n",xpad,ypad,c->fX[0],c->fY[0]);
1382 dx = xpad - c->fX[0];
1383 dy = ypad - c->fY[0];
1384 dr = TMath::Sqrt(dx*dx+dy*dy);
1389 fprintf(stderr," dr %f\n",dr);
1390 Int_t q=dig->Signal();
1391 if (dig->Physics() >= dig->Signal()) {
1392 c->fPhysicsMap[i]=2;
1393 } else if (dig->Physics() == 0) {
1394 c->fPhysicsMap[i]=0;
1395 } else c->fPhysicsMap[i]=1;
1396 c->fPeakSignal[cath]=q;
1397 c->fTracks[0]=dig->Hit();
1398 c->fTracks[1]=dig->Track(0);
1399 c->fTracks[2]=dig->Track(1);
1401 fprintf(stderr," c->fTracks[0] %d c->fTracks[1] %d\n",dig->Hit(),
1405 } // loop over digits
1407 // apply correction to the coordinate along the anode wire
1409 c->fX[cath]=fSeg[cath]->GetAnod(c->fX[cath]);
1412 void AliMUONClusterFinderVS::FindCluster(Int_t i, Int_t j, Int_t cath, AliMUONRawCluster &c){
1416 // Find a super cluster on both cathodes
1419 // Add i,j as element of the cluster
1422 Int_t idx = fHitMap[cath]->GetHitIndex(i,j);
1423 AliMUONDigit* dig = (AliMUONDigit*) fHitMap[cath]->GetHit(i,j);
1424 Int_t q=dig->Signal();
1425 Int_t theX=dig->PadX();
1426 Int_t theY=dig->PadY();
1428 if (q > TMath::Abs(c.fPeakSignal[0]) && q > TMath::Abs(c.fPeakSignal[1])) {
1429 c.fPeakSignal[cath]=q;
1430 c.fTracks[0]=dig->Hit();
1431 c.fTracks[1]=dig->Track(0);
1432 c.fTracks[2]=dig->Track(1);
1436 // Make sure that list of digits is ordered
1438 Int_t mu=c.fMultiplicity[cath];
1439 c.fIndexMap[mu][cath]=idx;
1441 if (dig->Physics() >= dig->Signal()) {
1442 c.fPhysicsMap[mu]=2;
1443 } else if (dig->Physics() == 0) {
1444 c.fPhysicsMap[mu]=0;
1445 } else c.fPhysicsMap[mu]=1;
1449 for (Int_t ind = mu-1; ind >= 0; ind--) {
1450 Int_t ist=(c.fIndexMap)[ind][cath];
1451 Int_t ql=fInput->Digit(cath, ist)->Signal();
1452 Int_t ix=fInput->Digit(cath, ist)->PadX();
1453 Int_t iy=fInput->Digit(cath, ist)->PadY();
1455 if (q>ql || (q==ql && theX > ix && theY < iy)) {
1456 c.fIndexMap[ind][cath]=idx;
1457 c.fIndexMap[ind+1][cath]=ist;
1465 c.fMultiplicity[cath]++;
1466 if (c.fMultiplicity[cath] >= 50 ) {
1467 printf("FindCluster - multiplicity >50 %d \n",c.fMultiplicity[0]);
1468 c.fMultiplicity[cath]=49;
1471 // Prepare center of gravity calculation
1473 fSeg[cath]->GetPadC(i, j, x, y, z);
1479 // Flag hit as "taken"
1480 fHitMap[cath]->FlagHit(i,j);
1482 // Now look recursively for all neighbours and pad hit on opposite cathode
1484 // Loop over neighbours
1488 Int_t xList[10], yList[10];
1489 fSeg[cath]->Neighbours(i,j,&nn,xList,yList);
1490 for (Int_t in=0; in<nn; in++) {
1494 if (fHitMap[cath]->TestHit(ix,iy)==kUnused) {
1496 printf("\n Neighbours %d %d %d", cath, ix, iy);
1497 FindCluster(ix, iy, cath, c);
1502 Int_t iXopp[50], iYopp[50];
1504 // Neighbours on opposite cathode
1505 // Take into account that several pads can overlap with the present pad
1506 Int_t isec=fSeg[cath]->Sector(i,j);
1512 dx = (fSeg[cath]->Dpx(isec))/2.;
1517 dy = (fSeg[cath]->Dpy(isec))/2;
1519 // loop over pad neighbours on opposite cathode
1520 for (fSeg[iop]->FirstPad(x, y, fZPlane, dx, dy);
1521 fSeg[iop]->MorePads();
1522 fSeg[iop]->NextPad())
1525 ix = fSeg[iop]->Ix(); iy = fSeg[iop]->Iy();
1526 if (fDebugLevel > 1)
1527 printf("\n ix, iy: %f %f %f %d %d %d", x,y,z,ix, iy, fSector);
1528 if (fHitMap[iop]->TestHit(ix,iy)==kUnused){
1531 if (fDebugLevel > 1)
1532 printf("\n Opposite %d %d %d", iop, ix, iy);
1535 } // Loop over pad neighbours
1536 // This had to go outside the loop since recursive calls inside the iterator are not possible
1539 for (jopp=0; jopp<nOpp; jopp++) {
1540 if (fHitMap[iop]->TestHit(iXopp[jopp],iYopp[jopp]) == kUnused)
1541 FindCluster(iXopp[jopp], iYopp[jopp], iop, c);
1545 //_____________________________________________________________________________
1547 void AliMUONClusterFinderVS::FindRawClusters()
1550 // MUON cluster finder from digits -- finds neighbours on both cathodes and
1551 // fills the tree with raw clusters
1554 // Return if no input datad available
1555 if (!fInput->NDigits(0) && !fInput->NDigits(1)) return;
1557 fSeg[0] = fInput->Segmentation(0);
1558 fSeg[1] = fInput->Segmentation(1);
1560 fHitMap[0] = new AliMUONHitMapA1(fSeg[0], fInput->Digits(0));
1561 fHitMap[1] = new AliMUONHitMapA1(fSeg[1], fInput->Digits(1));
1569 fHitMap[0]->FillHits();
1570 fHitMap[1]->FillHits();
1572 // Outer Loop over Cathodes
1573 for (cath=0; cath<2; cath++) {
1574 for (ndig=0; ndig<fInput->NDigits(cath); ndig++) {
1575 dig = fInput->Digit(cath, ndig);
1576 Int_t i=dig->PadX();
1577 Int_t j=dig->PadY();
1578 if (fHitMap[cath]->TestHit(i,j)==kUsed ||fHitMap[0]->TestHit(i,j)==kEmpty) {
1583 fprintf(stderr,"\n CATHODE %d CLUSTER %d\n",cath,ncls);
1584 AliMUONRawCluster c;
1585 c.fMultiplicity[0]=0;
1586 c.fMultiplicity[1]=0;
1587 c.fPeakSignal[cath]=dig->Signal();
1588 c.fTracks[0]=dig->Hit();
1589 c.fTracks[1]=dig->Track(0);
1590 c.fTracks[2]=dig->Track(1);
1591 // tag the beginning of cluster list in a raw cluster
1594 fSeg[cath]->GetPadC(i,j,xcu, ycu, fZPlane);
1595 fSector= fSeg[cath]->Sector(i,j)/100;
1597 printf("\n New Seed %d %d ", i,j);
1599 FindCluster(i,j,cath,c);
1600 // ^^^^^^^^^^^^^^^^^^^^^^^^
1601 // center of gravity
1604 c.fX[0]=fSeg[0]->GetAnod(c.fX[0]);
1608 c.fX[1]=fSeg[0]->GetAnod(c.fX[1]);
1615 fprintf(stderr,"\n Cathode 1 multiplicite %d X(CG) %f Y(CG) %f\n",
1616 c.fMultiplicity[0],c.fX[0],c.fY[0]);
1617 fprintf(stderr," Cathode 2 multiplicite %d X(CG) %f Y(CG) %f\n",
1618 c.fMultiplicity[1],c.fX[1],c.fY[1]);
1620 // Analyse cluster and decluster if necessary
1623 c.fNcluster[1]=fNRawClusters;
1624 c.fClusterType=c.PhysicsContribution();
1631 // reset Cluster object
1632 { // begin local scope
1633 for (int k=0;k<c.fMultiplicity[0];k++) c.fIndexMap[k][0]=0;
1634 } // end local scope
1636 { // begin local scope
1637 for (int k=0;k<c.fMultiplicity[1];k++) c.fIndexMap[k][1]=0;
1638 } // end local scope
1640 c.fMultiplicity[0]=c.fMultiplicity[0]=0;
1644 } // end loop cathodes
1649 Float_t AliMUONClusterFinderVS::SingleMathiesonFit(AliMUONRawCluster *c, Int_t cath)
1651 // Performs a single Mathieson fit on one cathode
1653 AliMUONClusterInput& clusterInput = *(AliMUONClusterInput::Instance());
1655 clusterInput.Fitter()->SetFCN(fcnS1);
1656 clusterInput.Fitter()->mninit(2,10,7);
1657 Double_t arglist[20];
1660 // Set starting values
1661 static Double_t vstart[2];
1666 // lower and upper limits
1667 static Double_t lower[2], upper[2];
1669 fSeg[cath]->GetPadI(c->fX[cath], c->fY[cath], fZPlane, ix, iy);
1670 Int_t isec=fSeg[cath]->Sector(ix, iy);
1671 lower[0]=vstart[0]-fSeg[cath]->Dpx(isec)/2;
1672 lower[1]=vstart[1]-fSeg[cath]->Dpy(isec)/2;
1674 upper[0]=lower[0]+fSeg[cath]->Dpx(isec);
1675 upper[1]=lower[1]+fSeg[cath]->Dpy(isec);
1678 static Double_t step[2]={0.0005, 0.0005};
1680 clusterInput.Fitter()->mnparm(0,"x1",vstart[0],step[0],lower[0],upper[0],ierflag);
1681 clusterInput.Fitter()->mnparm(1,"y1",vstart[1],step[1],lower[1],upper[1],ierflag);
1682 // ready for minimisation
1683 clusterInput.Fitter()->SetPrintLevel(-1+fDebugLevel);
1685 clusterInput.Fitter()->mnexcm("SET NOW", arglist, 0, ierflag);
1686 clusterInput.Fitter()->mnexcm("SET OUT", arglist, 0, ierflag);
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 Float_t AliMUONClusterFinderVS::CombiSingleMathiesonFit(AliMUONRawCluster *c)
1714 // Perform combined Mathieson fit on both cathode planes
1716 AliMUONClusterInput& clusterInput = *(AliMUONClusterInput::Instance());
1717 clusterInput.Fitter()->SetFCN(fcnCombiS1);
1718 clusterInput.Fitter()->mninit(2,10,7);
1719 Double_t arglist[20];
1722 static Double_t vstart[2];
1723 vstart[0]=fXInit[0];
1724 vstart[1]=fYInit[0];
1727 // lower and upper limits
1728 static Float_t lower[2], upper[2];
1730 fSeg[0]->GetPadI(fXInit[0], fYInit[0], fZPlane, ix, iy);
1731 isec=fSeg[0]->Sector(ix, iy);
1732 Float_t dpy=fSeg[0]->Dpy(isec);
1733 fSeg[1]->GetPadI(fXInit[0], fYInit[0], fZPlane, ix, iy);
1734 isec=fSeg[1]->Sector(ix, iy);
1735 Float_t dpx=fSeg[1]->Dpx(isec);
1738 Float_t xdum, ydum, zdum;
1740 // Find save upper and lower limits
1744 for (fSeg[1]->FirstPad(fXInit[0], fYInit[0], fZPlane, dpx, 0.);
1745 fSeg[1]->MorePads(); fSeg[1]->NextPad())
1747 ix=fSeg[1]->Ix(); iy=fSeg[1]->Iy();
1748 fSeg[1]->GetPadC(ix,iy, upper[0], ydum, zdum);
1749 if (icount ==0) lower[0]=upper[0];
1753 if (lower[0]>upper[0]) {xdum=lower[0]; lower[0]=upper[0]; upper[0]=xdum;}
1757 printf("\n single y %f %f", fXInit[0], fYInit[0]);
1759 for (fSeg[0]->FirstPad(fXInit[0], fYInit[0], fZPlane, 0., dpy);
1760 fSeg[0]->MorePads(); fSeg[0]->NextPad())
1762 ix=fSeg[0]->Ix(); iy=fSeg[0]->Iy();
1763 fSeg[0]->GetPadC(ix,iy,xdum,upper[1],zdum);
1764 if (icount ==0) lower[1]=upper[1];
1767 printf("\n upper lower %d %f %f", icount, upper[1], lower[1]);
1770 if (lower[1]>upper[1]) {xdum=lower[1]; lower[1]=upper[1]; upper[1]=xdum;}
1773 static Double_t step[2]={0.00001, 0.0001};
1775 clusterInput.Fitter()->mnparm(0,"x1",vstart[0],step[0],lower[0],upper[0],ierflag);
1776 clusterInput.Fitter()->mnparm(1,"y1",vstart[1],step[1],lower[1],upper[1],ierflag);
1777 // ready for minimisation
1778 clusterInput.Fitter()->SetPrintLevel(-1+fDebugLevel);
1780 clusterInput.Fitter()->mnexcm("SET NOW", arglist, 0, ierflag);
1781 clusterInput.Fitter()->mnexcm("SET OUT", arglist, 0, ierflag);
1785 clusterInput.Fitter()->mnexcm("SET NOGR", arglist, 0, ierflag);
1786 clusterInput.Fitter()->mnexcm("MIGRAD", arglist, 0, ierflag);
1787 clusterInput.Fitter()->mnexcm("EXIT" , arglist, 0, ierflag);
1788 Double_t fmin, fedm, errdef;
1789 Int_t npari, nparx, istat;
1791 clusterInput.Fitter()->mnstat(fmin, fedm, errdef, npari, nparx, istat);
1795 // Get fitted parameters
1796 Double_t xrec, yrec;
1798 Double_t epxz, b1, b2;
1800 clusterInput.Fitter()->mnpout(0, chname, xrec, epxz, b1, b2, ierflg);
1801 clusterInput.Fitter()->mnpout(1, chname, yrec, epxz, b1, b2, ierflg);
1807 Bool_t AliMUONClusterFinderVS::DoubleMathiesonFit(AliMUONRawCluster *c, Int_t cath)
1809 // Performs a double Mathieson fit on one cathode
1813 // Initialise global variables for fit
1814 AliMUONClusterInput& clusterInput = *(AliMUONClusterInput::Instance());
1815 clusterInput.Fitter()->SetFCN(fcnS2);
1816 clusterInput.Fitter()->mninit(5,10,7);
1817 Double_t arglist[20];
1820 // Set starting values
1821 static Double_t vstart[5];
1822 vstart[0]=fX[fIndLocal[0][cath]][cath];
1823 vstart[1]=fY[fIndLocal[0][cath]][cath];
1824 vstart[2]=fX[fIndLocal[1][cath]][cath];
1825 vstart[3]=fY[fIndLocal[1][cath]][cath];
1826 vstart[4]=Float_t(fQ[fIndLocal[0][cath]][cath])/
1827 Float_t(fQ[fIndLocal[0][cath]][cath]+fQ[fIndLocal[1][cath]][cath]);
1828 // lower and upper limits
1829 static Float_t lower[5], upper[5];
1830 Int_t isec=fSeg[cath]->Sector(fIx[fIndLocal[0][cath]][cath], fIy[fIndLocal[0][cath]][cath]);
1831 lower[0]=vstart[0]-fSeg[cath]->Dpx(isec);
1832 lower[1]=vstart[1]-fSeg[cath]->Dpy(isec);
1834 upper[0]=lower[0]+2.*fSeg[cath]->Dpx(isec);
1835 upper[1]=lower[1]+2.*fSeg[cath]->Dpy(isec);
1837 isec=fSeg[cath]->Sector(fIx[fIndLocal[1][cath]][cath], fIy[fIndLocal[1][cath]][cath]);
1838 lower[2]=vstart[2]-fSeg[cath]->Dpx(isec)/2;
1839 lower[3]=vstart[3]-fSeg[cath]->Dpy(isec)/2;
1841 upper[2]=lower[2]+fSeg[cath]->Dpx(isec);
1842 upper[3]=lower[3]+fSeg[cath]->Dpy(isec);
1847 static Double_t step[5]={0.0005, 0.0005, 0.0005, 0.0005, 0.0001};
1849 clusterInput.Fitter()->mnparm(0,"x1",vstart[0],step[0],lower[0],upper[0],ierflag);
1850 clusterInput.Fitter()->mnparm(1,"y1",vstart[1],step[1],lower[1],upper[1],ierflag);
1851 clusterInput.Fitter()->mnparm(2,"x2",vstart[2],step[2],lower[2],upper[2],ierflag);
1852 clusterInput.Fitter()->mnparm(3,"y2",vstart[3],step[3],lower[3],upper[3],ierflag);
1853 clusterInput.Fitter()->mnparm(4,"a0",vstart[4],step[4],lower[4],upper[4],ierflag);
1854 // ready for minimisation
1855 clusterInput.Fitter()->SetPrintLevel(-1+fDebugLevel);
1857 clusterInput.Fitter()->mnexcm("SET NOW", arglist, 0, ierflag);
1858 clusterInput.Fitter()->mnexcm("SET OUT", arglist, 0, ierflag);
1862 clusterInput.Fitter()->mnexcm("SET NOGR", arglist, 0, ierflag);
1863 clusterInput.Fitter()->mnexcm("MIGRAD", arglist, 0, ierflag);
1864 clusterInput.Fitter()->mnexcm("EXIT" , arglist, 0, ierflag);
1865 // Get fitted parameters
1866 Double_t xrec[2], yrec[2], qfrac;
1868 Double_t epxz, b1, b2;
1870 clusterInput.Fitter()->mnpout(0, chname, xrec[0], epxz, b1, b2, ierflg);
1871 clusterInput.Fitter()->mnpout(1, chname, yrec[0], epxz, b1, b2, ierflg);
1872 clusterInput.Fitter()->mnpout(2, chname, xrec[1], epxz, b1, b2, ierflg);
1873 clusterInput.Fitter()->mnpout(3, chname, yrec[1], epxz, b1, b2, ierflg);
1874 clusterInput.Fitter()->mnpout(4, chname, qfrac, epxz, b1, b2, ierflg);
1876 Double_t fmin, fedm, errdef;
1877 Int_t npari, nparx, istat;
1879 clusterInput.Fitter()->mnstat(fmin, fedm, errdef, npari, nparx, istat);
1884 Float_t AliMUONClusterFinderVS::CombiDoubleMathiesonFit(AliMUONRawCluster *c)
1887 // Perform combined double Mathieson fit on both cathode planes
1889 AliMUONClusterInput& clusterInput = *(AliMUONClusterInput::Instance());
1890 clusterInput.Fitter()->SetFCN(fcnCombiS2);
1891 clusterInput.Fitter()->mninit(6,10,7);
1892 Double_t arglist[20];
1895 // Set starting values
1896 static Double_t vstart[6];
1897 vstart[0]=fXInit[0];
1898 vstart[1]=fYInit[0];
1899 vstart[2]=fXInit[1];
1900 vstart[3]=fYInit[1];
1901 vstart[4]=fQrInit[0];
1902 vstart[5]=fQrInit[1];
1903 // lower and upper limits
1904 static Float_t lower[6], upper[6];
1908 fSeg[1]->GetPadI(fXInit[0], fYInit[0], fZPlane, ix, iy);
1909 isec=fSeg[1]->Sector(ix, iy);
1910 dpx=fSeg[1]->Dpx(isec);
1912 fSeg[0]->GetPadI(fXInit[0], fYInit[0], fZPlane, ix, iy);
1913 isec=fSeg[0]->Sector(ix, iy);
1914 dpy=fSeg[0]->Dpy(isec);
1918 Float_t xdum, ydum, zdum;
1920 printf("\n Cluster Finder: %f %f %f %f ", fXInit[0], fXInit[1],fYInit[0], fYInit[1] );
1922 // Find save upper and lower limits
1925 for (fSeg[1]->FirstPad(fXInit[0], fYInit[0], fZPlane, dpx, 0.);
1926 fSeg[1]->MorePads(); fSeg[1]->NextPad())
1928 ix=fSeg[1]->Ix(); iy=fSeg[1]->Iy();
1929 // if (fHitMap[1]->TestHit(ix, iy) == kEmpty) continue;
1930 fSeg[1]->GetPadC(ix,iy,upper[0],ydum,zdum);
1931 if (icount ==0) lower[0]=upper[0];
1934 if (lower[0]>upper[0]) {xdum=lower[0]; lower[0]=upper[0]; upper[0]=xdum;}
1935 // vstart[0] = 0.5*(lower[0]+upper[0]);
1940 for (fSeg[0]->FirstPad(fXInit[0], fYInit[0], fZPlane, 0., dpy);
1941 fSeg[0]->MorePads(); fSeg[0]->NextPad())
1943 ix=fSeg[0]->Ix(); iy=fSeg[0]->Iy();
1944 // if (fHitMap[0]->TestHit(ix, iy) == kEmpty) continue;
1945 fSeg[0]->GetPadC(ix,iy,xdum,upper[1],zdum);
1946 if (icount ==0) lower[1]=upper[1];
1950 if (lower[1]>upper[1]) {xdum=lower[1]; lower[1]=upper[1]; upper[1]=xdum;}
1951 // vstart[1] = 0.5*(lower[1]+upper[1]);
1954 fSeg[1]->GetPadI(fXInit[1], fYInit[1], fZPlane, ix, iy);
1955 isec=fSeg[1]->Sector(ix, iy);
1956 dpx=fSeg[1]->Dpx(isec);
1957 fSeg[0]->GetPadI(fXInit[1], fYInit[1], fZPlane, ix, iy);
1958 isec=fSeg[0]->Sector(ix, iy);
1959 dpy=fSeg[0]->Dpy(isec);
1962 // Find save upper and lower limits
1966 for (fSeg[1]->FirstPad(fXInit[1], fYInit[1], fZPlane, dpx, 0);
1967 fSeg[1]->MorePads(); fSeg[1]->NextPad())
1969 ix=fSeg[1]->Ix(); iy=fSeg[1]->Iy();
1970 // if (fHitMap[1]->TestHit(ix, iy) == kEmpty) continue;
1971 fSeg[1]->GetPadC(ix,iy,upper[2],ydum,zdum);
1972 if (icount ==0) lower[2]=upper[2];
1975 if (lower[2]>upper[2]) {xdum=lower[2]; lower[2]=upper[2]; upper[2]=xdum;}
1976 // vstart[2] = 0.5*(lower[2]+upper[2]);
1980 for (fSeg[0]->FirstPad(fXInit[1], fYInit[1], fZPlane, 0, dpy);
1981 fSeg[0]-> MorePads(); fSeg[0]->NextPad())
1983 ix=fSeg[0]->Ix(); iy=fSeg[0]->Iy();
1984 // if (fHitMap[0]->TestHit(ix, iy) != kEmpty) continue;
1986 fSeg[0]->GetPadC(ix,iy,xdum,upper[3],zdum);
1987 if (icount ==0) lower[3]=upper[3];
1991 if (lower[3]>upper[3]) {xdum=lower[3]; lower[3]=upper[3]; upper[3]=xdum;}
1993 // vstart[3] = 0.5*(lower[3]+upper[3]);
2001 static Double_t step[6]={0.0005, 0.0005, 0.0005, 0.0005, 0.001, 0.001};
2002 clusterInput.Fitter()->mnparm(0,"x1",vstart[0],step[0],lower[0],upper[0],ierflag);
2003 clusterInput.Fitter()->mnparm(1,"y1",vstart[1],step[1],lower[1],upper[1],ierflag);
2004 clusterInput.Fitter()->mnparm(2,"x2",vstart[2],step[2],lower[2],upper[2],ierflag);
2005 clusterInput.Fitter()->mnparm(3,"y2",vstart[3],step[3],lower[3],upper[3],ierflag);
2006 clusterInput.Fitter()->mnparm(4,"a0",vstart[4],step[4],lower[4],upper[4],ierflag);
2007 clusterInput.Fitter()->mnparm(5,"a1",vstart[5],step[5],lower[5],upper[5],ierflag);
2008 // ready for minimisation
2009 clusterInput.Fitter()->SetPrintLevel(-1+fDebugLevel);
2011 clusterInput.Fitter()->mnexcm("SET NOW", arglist, 0, ierflag);
2012 clusterInput.Fitter()->mnexcm("SET OUT", arglist, 0, ierflag);
2016 clusterInput.Fitter()->mnexcm("SET NOGR", arglist, 0, ierflag);
2017 clusterInput.Fitter()->mnexcm("MIGRAD", arglist, 0, ierflag);
2018 clusterInput.Fitter()->mnexcm("EXIT" , arglist, 0, ierflag);
2019 // Get fitted parameters
2021 Double_t epxz, b1, b2;
2023 clusterInput.Fitter()->mnpout(0, chname, fXFit[0], epxz, b1, b2, ierflg);
2024 clusterInput.Fitter()->mnpout(1, chname, fYFit[0], epxz, b1, b2, ierflg);
2025 clusterInput.Fitter()->mnpout(2, chname, fXFit[1], epxz, b1, b2, ierflg);
2026 clusterInput.Fitter()->mnpout(3, chname, fYFit[1], epxz, b1, b2, ierflg);
2027 clusterInput.Fitter()->mnpout(4, chname, fQrFit[0], epxz, b1, b2, ierflg);
2028 clusterInput.Fitter()->mnpout(5, chname, fQrFit[1], epxz, b1, b2, ierflg);
2030 Double_t fmin, fedm, errdef;
2031 Int_t npari, nparx, istat;
2033 clusterInput.Fitter()->mnstat(fmin, fedm, errdef, npari, nparx, istat);
2041 void AliMUONClusterFinderVS::Split(AliMUONRawCluster* c)
2044 // One cluster for each maximum
2047 AliMUONClusterInput& clusterInput = *(AliMUONClusterInput::Instance());
2048 for (j=0; j<2; j++) {
2049 AliMUONRawCluster cnew;
2050 cnew.fGhost=c->fGhost;
2051 for (cath=0; cath<2; cath++) {
2052 cnew.fChi2[cath]=fChi2[0];
2053 // ?? why not cnew.fChi2[cath]=fChi2[cath];
2056 cnew.fNcluster[0]=-1;
2057 cnew.fNcluster[1]=fNRawClusters;
2059 cnew.fNcluster[0]=fNPeaks;
2060 cnew.fNcluster[1]=0;
2062 cnew.fMultiplicity[cath]=0;
2063 cnew.fX[cath]=Float_t(fXFit[j]);
2064 cnew.fY[cath]=Float_t(fYFit[j]);
2065 cnew.fZ[cath]=fZPlane;
2067 cnew.fQ[cath]=Int_t(clusterInput.TotalCharge(cath)*fQrFit[cath]);
2069 cnew.fQ[cath]=Int_t(clusterInput.TotalCharge(cath)*(1-fQrFit[cath]));
2071 fSeg[cath]->SetHit(fXFit[j],fYFit[j],fZPlane);
2072 for (i=0; i<fMul[cath]; i++) {
2073 cnew.fIndexMap[cnew.fMultiplicity[cath]][cath]=
2074 c->fIndexMap[i][cath];
2075 fSeg[cath]->SetPad(fIx[i][cath], fIy[i][cath]);
2076 Float_t q1=fInput->Response()->IntXY(fSeg[cath]);
2077 cnew.fContMap[i][cath]
2078 =(q1*Float_t(cnew.fQ[cath]))/Float_t(fQ[i][cath]);
2079 cnew.fMultiplicity[cath]++;
2081 FillCluster(&cnew,0,cath);
2084 cnew.fClusterType=cnew.PhysicsContribution();
2085 if (cnew.fQ[0]>0 && cnew.fQ[1]>0) AddRawCluster(cnew);
2092 // Minimisation functions
2094 void fcnS1(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag)
2096 AliMUONClusterInput& clusterInput = *(AliMUONClusterInput::Instance());
2103 for (i=0; i<clusterInput.Nmul(0); i++) {
2104 Float_t q0=clusterInput.Charge(i,0);
2105 Float_t q1=clusterInput.DiscrChargeS1(i,par);
2114 void fcnCombiS1(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag)
2116 AliMUONClusterInput& clusterInput = *(AliMUONClusterInput::Instance());
2123 for (cath=0; cath<2; cath++) {
2124 for (i=0; i<clusterInput.Nmul(cath); i++) {
2125 Float_t q0=clusterInput.Charge(i,cath);
2126 Float_t q1=clusterInput.DiscrChargeCombiS1(i,par,cath);
2137 void fcnS2(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag)
2139 AliMUONClusterInput& clusterInput = *(AliMUONClusterInput::Instance());
2146 for (i=0; i<clusterInput.Nmul(0); i++) {
2148 Float_t q0=clusterInput.Charge(i,0);
2149 Float_t q1=clusterInput.DiscrChargeS2(i,par);
2159 void fcnCombiS2(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag)
2161 AliMUONClusterInput& clusterInput = *(AliMUONClusterInput::Instance());
2167 for (cath=0; cath<2; cath++) {
2168 for (i=0; i<clusterInput.Nmul(cath); i++) {
2169 Float_t q0=clusterInput.Charge(i,cath);
2170 Float_t q1=clusterInput.DiscrChargeCombiS2(i,par,cath);
2180 void AliMUONClusterFinderVS::AddRawCluster(const AliMUONRawCluster c)
2183 // Add a raw cluster copy to the list
2185 AliMUON *pMUON=(AliMUON*)gAlice->GetModule("MUON");
2186 pMUON->AddRawCluster(fInput->Chamber(),c);
2189 fprintf(stderr,"\nfNRawClusters %d\n",fNRawClusters);
2192 Bool_t AliMUONClusterFinderVS::TestTrack(Int_t t) {
2193 // Test if track was user selected
2194 if (fTrack[0]==-1 || fTrack[1]==-1) {
2196 } else if (t==fTrack[0] || t==fTrack[1]) {
2203 AliMUONClusterFinderVS& AliMUONClusterFinderVS
2204 ::operator = (const AliMUONClusterFinderVS& rhs)
2206 // Dummy assignment operator