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.20 2001/04/12 12:22:26 morsch
18 - some numerical problems caused by pad staggering cured.
19 - treatment of 1-2 and 2-1 ghosts
20 - debuglevel > 1 prints introduced
22 Revision 1.19 2001/03/20 13:32:10 egangler
23 Code introduced to remove ghosts with the charge correlation between the 2
24 cathods. A chi2 is performed for the 2 possibilities.
25 If one gets good chi2 (with respect to the fGhostChi2Cut parameter) and the
26 other wrong chi2, the ambiguity is solved
27 If both gets good or both bad chi2, no choice is made
28 By default the fGhostChi2Cut parameter is set to solve about 70% of ghost
29 problems with about 2% errors, with the current version of the code.
32 fGhostChi2Cut is in AliMUONClusterFinderVS, with setters and getters.
33 a fDebugLevel was also introduced to switch off some of the output.
34 When an ambiguity is detected and not solved, the fGhost word in
35 AliMUONRawCluster is set to 1 or 2, depending whether both charge chi2 are
37 a DumpIndex method was also added in AliMUONRawCluster to produce a printout
41 By default, the code makes ghost check. If you want previous behaviour,
42 put in MUONrawclusters the value of SetGhostChi2Cut to infinity (1e.6) is
45 Revision 1.18 2001/01/26 21:37:53 morsch
46 Use access functions to AliMUONDigit member data.
48 Revision 1.17 2001/01/23 18:58:19 hristov
49 Initialisation of some pointers
51 Revision 1.16 2000/12/21 23:27:30 morsch
52 Error in argument list of AddRawCluster corrected.
54 Revision 1.15 2000/12/21 22:14:38 morsch
55 Clean-up of coding rule violations.
57 Revision 1.14 2000/10/23 16:03:45 morsch
58 Correct z-position of all clusters created "on the flight".
60 Revision 1.13 2000/10/23 13:38:23 morsch
61 Set correct z-coordinate when cluster is split.
63 Revision 1.12 2000/10/18 11:42:06 morsch
64 - AliMUONRawCluster contains z-position.
65 - Some clean-up of useless print statements during initialisations.
67 Revision 1.11 2000/10/06 09:04:05 morsch
68 - Dummy z-arguments in GetPadI, SetHit, FirstPad replaced by real z-coordinate
69 to make code work with slat chambers (AM)
70 - Replace GetPadI calls with unchecked x,y coordinates by pad iterator calls wherever possible.
72 Revision 1.10 2000/10/03 13:51:57 egangler
73 Removal of useless dependencies via forward declarations
75 Revision 1.9 2000/10/02 16:58:29 egangler
76 Cleaning of the code :
79 -> some useless includes removed or replaced by "class" statement
81 Revision 1.8 2000/07/03 11:54:57 morsch
82 AliMUONSegmentation and AliMUONHitMap have been replaced by AliSegmentation and AliHitMap in STEER
83 The methods GetPadIxy and GetPadXxy of AliMUONSegmentation have changed name to GetPadI and GetPadC.
85 Revision 1.7 2000/06/28 15:16:35 morsch
86 (1) Client code adapted to new method signatures in AliMUONSegmentation (see comments there)
87 to allow development of slat-muon chamber simulation and reconstruction code in the MUON
88 framework. The changes should have no side effects (mostly dummy arguments).
89 (2) Hit disintegration uses 3-dim hit coordinates to allow simulation
90 of chambers with overlapping modules (MakePadHits, Disintegration).
92 Revision 1.6 2000/06/28 12:19:18 morsch
93 More consequent seperation of global input data services (AliMUONClusterInput singleton) and the
94 cluster and hit reconstruction algorithms in AliMUONClusterFinderVS.
95 AliMUONClusterFinderVS becomes the base class for clustering and hit reconstruction.
96 It requires two cathode planes. Small modifications in the code will make it usable for
97 one cathode plane and, hence, more general (for test beam data).
98 AliMUONClusterFinder is now obsolete.
100 Revision 1.5 2000/06/28 08:06:10 morsch
101 Avoid global variables in AliMUONClusterFinderVS by seperating the input data for the fit from the
102 algorithmic part of the class. Input data resides inside the AliMUONClusterInput singleton.
103 It also naturally takes care of the TMinuit instance.
105 Revision 1.4 2000/06/27 16:18:47 gosset
106 Finally correct implementation of xm, ym, ixm, iym sizes
107 when at least three local maxima on cathode 1 or on cathode 2
109 Revision 1.3 2000/06/22 14:02:45 morsch
110 Parameterised size of xm[], ym[], ixm[], iym[] correctly implemented (PH)
111 Some HP scope problems corrected (PH)
113 Revision 1.2 2000/06/15 07:58:48 morsch
114 Code from MUON-dev joined
116 Revision 1.1.2.3 2000/06/09 21:58:33 morsch
117 Most coding rule violations corrected.
119 Revision 1.1.2.2 2000/02/15 08:33:52 morsch
120 Error in calculation of contribution map for double clusters (Split method) corrected (A.M.)
121 Error in determination of track list for double cluster (FillCluster method) corrected (A.M.)
122 Revised and extended SplitByLocalMaxima method (Isabelle Chevrot):
123 - For clusters with more than 2 maxima on one of the cathode planes all valid
124 combinations of maxima on the two cathodes are preserved. The position of the maxima is
125 taken as the hit position.
126 - New FillCluster method with 2 arguments to find tracks associated to the clusters
127 defined above added. (Method destinction by argument list not very elegant in this case,
128 should be revides (A.M.)
129 - Bug in if-statement to handle maximum 1 maximum per plane corrected
130 - Two cluster per cathode but only 1 combination valid is handled.
131 - More rigerous treatment of 1-2 and 2-1 combinations of maxima.
135 #include "AliMUONClusterFinderVS.h"
136 #include "AliMUONDigit.h"
137 #include "AliMUONRawCluster.h"
138 #include "AliSegmentation.h"
139 #include "AliMUONResponse.h"
140 #include "AliMUONClusterInput.h"
141 #include "AliMUONHitMapA1.h"
150 #include <TPostScript.h>
155 #include <iostream.h>
157 //_____________________________________________________________________
158 // This function is minimized in the double-Mathieson fit
159 void fcnS2(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag);
160 void fcnS1(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag);
161 void fcnCombiS1(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag);
162 void fcnCombiS2(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag);
164 ClassImp(AliMUONClusterFinderVS)
166 AliMUONClusterFinderVS::AliMUONClusterFinderVS()
168 // Default constructor
169 fInput=AliMUONClusterInput::Instance();
172 fTrack[0]=fTrack[1]=-1;
173 fDebugLevel = 0; // make silent default
174 fGhostChi2Cut = 1e6; // nothing done by default
177 for(Int_t i=0; i<100; i++) {
178 for (Int_t j=0; j<2; j++) {
184 AliMUONClusterFinderVS::AliMUONClusterFinderVS(
185 const AliMUONClusterFinderVS & clusterFinder)
187 // Dummy copy Constructor
191 void AliMUONClusterFinderVS::Decluster(AliMUONRawCluster *cluster)
193 // Decluster by local maxima
194 SplitByLocalMaxima(cluster);
197 void AliMUONClusterFinderVS::SplitByLocalMaxima(AliMUONRawCluster *c)
199 // Split complex cluster by local maxima
202 fInput->SetCluster(c);
204 fMul[0]=c->fMultiplicity[0];
205 fMul[1]=c->fMultiplicity[1];
208 // dump digit information into arrays
213 for (cath=0; cath<2; cath++) {
215 for (i=0; i<fMul[cath]; i++)
218 fDig[i][cath]=fInput->Digit(cath, c->fIndexMap[i][cath]);
220 fIx[i][cath]= fDig[i][cath]->PadX();
221 fIy[i][cath]= fDig[i][cath]->PadY();
223 fQ[i][cath] = fDig[i][cath]->Signal();
224 // pad centre coordinates
226 GetPadC(fIx[i][cath], fIy[i][cath], fX[i][cath], fY[i][cath], fZ[i][cath]);
227 } // loop over cluster digits
228 } // loop over cathodes
234 // Initialise and perform mathieson fits
235 Float_t chi2, oldchi2;
236 // ++++++++++++++++++*************+++++++++++++++++++++
237 // (1) No more than one local maximum per cathode plane
238 // +++++++++++++++++++++++++++++++*************++++++++
239 if ((fNLocal[0]==1 && (fNLocal[1]==0 || fNLocal[1]==1)) ||
240 (fNLocal[0]==0 && fNLocal[1]==1)) {
241 // Perform combined single Mathieson fit
242 // Initial values for coordinates (x,y)
244 // One local maximum on cathodes 1 and 2 (X->cathode 2, Y->cathode 1)
245 if (fNLocal[0]==1 && fNLocal[1]==1) {
248 // One local maximum on cathode 1 (X,Y->cathode 1)
249 } else if (fNLocal[0]==1) {
252 // One local maximum on cathode 2 (X,Y->cathode 2)
258 fprintf(stderr,"\n cas (1) CombiSingleMathiesonFit(c)\n");
259 chi2=CombiSingleMathiesonFit(c);
260 // Int_t ndf = fgNbins[0]+fgNbins[1]-2;
261 // Float_t prob = TMath::Prob(Double_t(chi2),ndf);
262 // prob1->Fill(prob);
263 // chi2_1->Fill(chi2);
266 fprintf(stderr," chi2 %f ",chi2);
276 c->fX[0]=fSeg[0]->GetAnod(c->fX[0]);
277 c->fX[1]=fSeg[1]->GetAnod(c->fX[1]);
279 // If reasonable chi^2 add result to the list of rawclusters
282 // If not try combined double Mathieson Fit
284 fprintf(stderr," MAUVAIS CHI2 !!!\n");
285 if (fNLocal[0]==1 && fNLocal[1]==1) {
286 fXInit[0]=fX[fIndLocal[0][1]][1];
287 fYInit[0]=fY[fIndLocal[0][0]][0];
288 fXInit[1]=fX[fIndLocal[0][1]][1];
289 fYInit[1]=fY[fIndLocal[0][0]][0];
290 } else if (fNLocal[0]==1) {
291 fXInit[0]=fX[fIndLocal[0][0]][0];
292 fYInit[0]=fY[fIndLocal[0][0]][0];
293 fXInit[1]=fX[fIndLocal[0][0]][0];
294 fYInit[1]=fY[fIndLocal[0][0]][0];
296 fXInit[0]=fX[fIndLocal[0][1]][1];
297 fYInit[0]=fY[fIndLocal[0][1]][1];
298 fXInit[1]=fX[fIndLocal[0][1]][1];
299 fYInit[1]=fY[fIndLocal[0][1]][1];
302 // Initial value for charge ratios
306 fprintf(stderr,"\n cas (1) CombiDoubleMathiesonFit(c)\n");
307 chi2=CombiDoubleMathiesonFit(c);
308 // Int_t ndf = fgNbins[0]+fgNbins[1]-6;
309 // Float_t prob = TMath::Prob(chi2,ndf);
310 // prob2->Fill(prob);
311 // chi2_2->Fill(chi2);
313 // Was this any better ??
314 fprintf(stderr," Old and new chi2 %f %f ", oldchi2, chi2);
315 if (fFitStat!=0 && chi2>0 && (2.*chi2 < oldchi2)) {
316 fprintf(stderr," Split\n");
317 // Split cluster into two according to fit result
320 fprintf(stderr," Don't Split\n");
326 // +++++++++++++++++++++++++++++++++++++++
327 // (2) Two local maxima per cathode plane
328 // +++++++++++++++++++++++++++++++++++++++
329 } else if (fNLocal[0]==2 && fNLocal[1]==2) {
331 // Let's look for ghosts first
333 Float_t xm[4][2], ym[4][2];
334 Float_t dpx, dpy, dx, dy;
335 Int_t ixm[4][2], iym[4][2];
336 Int_t isec, im1, im2, ico;
338 // Form the 2x2 combinations
339 // 0-0, 0-1, 1-0, 1-1
341 for (im1=0; im1<2; im1++) {
342 for (im2=0; im2<2; im2++) {
343 xm[ico][0]=fX[fIndLocal[im1][0]][0];
344 ym[ico][0]=fY[fIndLocal[im1][0]][0];
345 xm[ico][1]=fX[fIndLocal[im2][1]][1];
346 ym[ico][1]=fY[fIndLocal[im2][1]][1];
348 ixm[ico][0]=fIx[fIndLocal[im1][0]][0];
349 iym[ico][0]=fIy[fIndLocal[im1][0]][0];
350 ixm[ico][1]=fIx[fIndLocal[im2][1]][1];
351 iym[ico][1]=fIy[fIndLocal[im2][1]][1];
355 // ico = 0 : first local maximum on cathodes 1 and 2
356 // ico = 1 : fisrt local maximum on cathode 1 and second on cathode 2
357 // ico = 2 : second local maximum on cathode 1 and first on cathode 1
358 // ico = 3 : second local maximum on cathodes 1 and 2
360 // Analyse the combinations and keep those that are possible !
361 // For each combination check consistency in x and y
364 Float_t dr[4] = {1.e4, 1.e4, 1.e4, 1.e4};
367 // In case of staggering maxima are displaced by exactly half the pad-size in y.
368 // We have to take into account the numerical precision in the consistency check;
371 for (ico=0; ico<4; ico++) {
372 accepted[ico]=kFALSE;
373 // cathode one: x-coordinate
374 isec=fSeg[0]->Sector(ixm[ico][0], iym[ico][0]);
375 dpx=fSeg[0]->Dpx(isec)/2.;
376 dx=TMath::Abs(xm[ico][0]-xm[ico][1]);
377 // cathode two: y-coordinate
378 isec=fSeg[1]->Sector(ixm[ico][1], iym[ico][1]);
379 dpy=fSeg[1]->Dpy(isec)/2.;
380 dy=TMath::Abs(ym[ico][0]-ym[ico][1]);
382 printf("\n %i %f %f %f %f %f %f \n", ico, ym[ico][0], ym[ico][1], dy, dpy, dx, dpx );
383 if ((dx <= dpx) && (dy <= dpy+eps)) {
386 dr[ico] = TMath::Sqrt(dx*dx+dy*dy);
390 accepted[ico]=kFALSE;
393 printf("\n iacc= %d:\n", iacc);
395 if (accepted[0] && accepted[1]) {
396 if (dr[0] >= dr[1]) {
403 if (accepted[2] && accepted[3]) {
404 if (dr[2] >= dr[3]) {
411 // eliminate one candidate
415 for (ico=0; ico<4; ico++) {
416 if (accepted[ico] && dr[ico] > drmax) {
422 accepted[icobad] = kFALSE;
428 printf("\n iacc= %d:\n", iacc);
431 fprintf(stderr,"\n iacc=2: No problem ! \n");
432 } else if (iacc==4) {
433 fprintf(stderr,"\n iacc=4: Ok, but ghost problem !!! \n");
434 } else if (iacc==0) {
435 fprintf(stderr,"\n iacc=0: I don't know what to do with this !!!!!!!!! \n");
439 // Initial value for charge ratios
440 fQrInit[0]=Float_t(fQ[fIndLocal[0][0]][0])/
441 Float_t(fQ[fIndLocal[0][0]][0]+fQ[fIndLocal[1][0]][0]);
442 fQrInit[1]=Float_t(fQ[fIndLocal[0][1]][1])/
443 Float_t(fQ[fIndLocal[0][1]][1]+fQ[fIndLocal[1][1]][1]);
445 // ******* iacc = 0 *******
446 // No combinations found between the 2 cathodes
447 // We keep the center of gravity of the cluster
452 // ******* iacc = 1 *******
453 // Only one combination found between the 2 cathodes
455 // Initial values for the 2 maxima (x,y)
457 // 1 maximum is initialised with the maximum of the combination found (X->cathode 2, Y->cathode 1)
458 // 1 maximum is initialised with the other maximum of the first cathode
460 fprintf(stderr,"ico=0\n");
465 } else if (accepted[1]){
466 fprintf(stderr,"ico=1\n");
471 } else if (accepted[2]){
472 fprintf(stderr,"ico=2\n");
477 } else if (accepted[3]){
478 fprintf(stderr,"ico=3\n");
485 fprintf(stderr,"\n cas (2) CombiDoubleMathiesonFit(c)\n");
486 chi2=CombiDoubleMathiesonFit(c);
487 // Int_t ndf = fgNbins[0]+fgNbins[1]-6;
488 // Float_t prob = TMath::Prob(chi2,ndf);
489 // prob2->Fill(prob);
490 // chi2_2->Fill(chi2);
492 fprintf(stderr," chi2 %f\n",chi2);
494 // If reasonable chi^2 add result to the list of rawclusters
499 // 1 maximum is initialised with the maximum of the combination found (X->cathode 2, Y->cathode 1)
500 // 1 maximum is initialised with the other maximum of the second cathode
502 fprintf(stderr,"ico=0\n");
507 } else if (accepted[1]){
508 fprintf(stderr,"ico=1\n");
513 } else if (accepted[2]){
514 fprintf(stderr,"ico=2\n");
519 } else if (accepted[3]){
520 fprintf(stderr,"ico=3\n");
527 fprintf(stderr,"\n cas (2) CombiDoubleMathiesonFit(c)\n");
528 chi2=CombiDoubleMathiesonFit(c);
529 // Int_t ndf = fgNbins[0]+fgNbins[1]-6;
530 // Float_t prob = TMath::Prob(chi2,ndf);
531 // prob2->Fill(prob);
532 // chi2_2->Fill(chi2);
534 fprintf(stderr," chi2 %f\n",chi2);
536 // If reasonable chi^2 add result to the list of rawclusters
540 //We keep only the combination found (X->cathode 2, Y->cathode 1)
541 for (Int_t ico=0; ico<2; ico++) {
543 AliMUONRawCluster cnew;
545 for (cath=0; cath<2; cath++) {
546 cnew.fX[cath]=Float_t(xm[ico][1]);
547 cnew.fY[cath]=Float_t(ym[ico][0]);
548 cnew.fZ[cath]=fZPlane;
550 cnew.fMultiplicity[cath]=c->fMultiplicity[cath];
551 for (i=0; i<fMul[cath]; i++) {
552 cnew.fIndexMap[i][cath]=c->fIndexMap[i][cath];
553 fSeg[cath]->SetPad(fIx[i][cath], fIy[i][cath]);
555 fprintf(stderr,"\nRawCluster %d cath %d\n",ico,cath);
556 fprintf(stderr,"mult_av %d\n",c->fMultiplicity[cath]);
557 FillCluster(&cnew,cath);
559 cnew.fClusterType=cnew.PhysicsContribution();
568 // ******* iacc = 2 *******
569 // Two combinations found between the 2 cathodes
571 // Was the same maximum taken twice
572 if ((accepted[0]&&accepted[1]) || (accepted[2]&&accepted[3])) {
573 fprintf(stderr,"\n Maximum taken twice !!!\n");
575 // Have a try !! with that
576 if (accepted[0]&&accepted[3]) {
588 fprintf(stderr,"\n cas (2) CombiDoubleMathiesonFit(c)\n");
589 chi2=CombiDoubleMathiesonFit(c);
590 // Int_t ndf = fgNbins[0]+fgNbins[1]-6;
591 // Float_t prob = TMath::Prob(chi2,ndf);
592 // prob2->Fill(prob);
593 // chi2_2->Fill(chi2);
597 // No ghosts ! No Problems ! - Perform one fit only !
598 if (accepted[0]&&accepted[3]) {
610 fprintf(stderr,"\n cas (2) CombiDoubleMathiesonFit(c)\n");
611 chi2=CombiDoubleMathiesonFit(c);
612 // Int_t ndf = fgNbins[0]+fgNbins[1]-6;
613 // Float_t prob = TMath::Prob(chi2,ndf);
614 // prob2->Fill(prob);
615 // chi2_2->Fill(chi2);
617 fprintf(stderr," chi2 %f\n",chi2);
621 // ******* iacc = 4 *******
622 // Four combinations found between the 2 cathodes
624 } else if (iacc==4) {
625 // Perform fits for the two possibilities !!
626 // Accept if charges are compatible on both cathodes
627 // If none are compatible, keep everything
633 fprintf(stderr,"\n cas (2) CombiDoubleMathiesonFit(c)\n");
634 chi2=CombiDoubleMathiesonFit(c);
635 // Int_t ndf = fgNbins[0]+fgNbins[1]-6;
636 // Float_t prob = TMath::Prob(chi2,ndf);
637 // prob2->Fill(prob);
638 // chi2_2->Fill(chi2);
640 fprintf(stderr," chi2 %f\n",chi2);
641 // store results of fit and postpone decision
642 Double_t sXFit[2],sYFit[2],sQrFit[2];
644 for (Int_t i=0;i<2;i++) {
655 fprintf(stderr,"\n cas (2) CombiDoubleMathiesonFit(c)\n");
656 chi2=CombiDoubleMathiesonFit(c);
657 // ndf = fgNbins[0]+fgNbins[1]-6;
658 // prob = TMath::Prob(chi2,ndf);
659 // prob2->Fill(prob);
660 // chi2_2->Fill(chi2);
662 fprintf(stderr," chi2 %f\n",chi2);
663 // We have all informations to perform the decision
664 // Compute the chi2 for the 2 possibilities
665 Float_t chi2fi,chi2si,chi2f,chi2s;
667 chi2f = (TMath::Log(fInput->TotalCharge(0)*fQrFit[0]
668 / (fInput->TotalCharge(1)*fQrFit[1]) )
669 / fInput->Response()->ChargeCorrel() );
671 chi2fi = (TMath::Log(fInput->TotalCharge(0)*(1-fQrFit[0])
672 / (fInput->TotalCharge(1)*(1-fQrFit[1])) )
673 / fInput->Response()->ChargeCorrel() );
674 chi2f += chi2fi*chi2fi;
676 chi2s = (TMath::Log(fInput->TotalCharge(0)*sQrFit[0]
677 / (fInput->TotalCharge(1)*sQrFit[1]) )
678 / fInput->Response()->ChargeCorrel() );
680 chi2si = (TMath::Log(fInput->TotalCharge(0)*(1-sQrFit[0])
681 / (fInput->TotalCharge(1)*(1-sQrFit[1])) )
682 / fInput->Response()->ChargeCorrel() );
683 chi2s += chi2si*chi2si;
685 // usefull to store the charge matching chi2 in the cluster
686 // fChi2[0]=sChi2[1]=chi2f;
687 // fChi2[1]=sChi2[0]=chi2s;
689 if (chi2f<=fGhostChi2Cut && chi2s<=fGhostChi2Cut)
691 if (chi2f>fGhostChi2Cut && chi2s>fGhostChi2Cut) {
697 if (chi2f<=fGhostChi2Cut)
699 if (chi2s<=fGhostChi2Cut) {
700 // retreive saved values
701 for (Int_t i=0;i<2;i++) {
712 } else if (fNLocal[0]==2 && fNLocal[1]==1) {
713 // ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
714 // (3) Two local maxima on cathode 1 and one maximum on cathode 2
715 // ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
717 Float_t xm[4][2], ym[4][2];
718 Float_t dpx, dpy, dx, dy;
719 Int_t ixm[4][2], iym[4][2];
720 Int_t isec, im1, ico;
722 // Form the 2x2 combinations
723 // 0-0, 0-1, 1-0, 1-1
725 for (im1=0; im1<2; im1++) {
726 xm[ico][0]=fX[fIndLocal[im1][0]][0];
727 ym[ico][0]=fY[fIndLocal[im1][0]][0];
728 xm[ico][1]=fX[fIndLocal[0][1]][1];
729 ym[ico][1]=fY[fIndLocal[0][1]][1];
731 ixm[ico][0]=fIx[fIndLocal[im1][0]][0];
732 iym[ico][0]=fIy[fIndLocal[im1][0]][0];
733 ixm[ico][1]=fIx[fIndLocal[0][1]][1];
734 iym[ico][1]=fIy[fIndLocal[0][1]][1];
737 // ico = 0 : first local maximum on cathodes 1 and 2
738 // ico = 1 : second local maximum on cathode 1 and first on cathode 2
740 // Analyse the combinations and keep those that are possible !
741 // For each combination check consistency in x and y
745 // In case of staggering maxima are displaced by exactly half the pad-size in y.
746 // We have to take into account the numerical precision in the consistency check;
750 for (ico=0; ico<2; ico++) {
751 accepted[ico]=kFALSE;
752 isec=fSeg[0]->Sector(ixm[ico][0], iym[ico][0]);
753 dpx=fSeg[0]->Dpx(isec)/2.;
754 dx=TMath::Abs(xm[ico][0]-xm[ico][1]);
755 isec=fSeg[1]->Sector(ixm[ico][1], iym[ico][1]);
756 dpy=fSeg[1]->Dpy(isec)/2.;
757 dy=TMath::Abs(ym[ico][0]-ym[ico][1]);
759 printf("\n %i %f %f %f %f \n", ico, ym[ico][0], ym[ico][1], dy, dpy );
760 if ((dx <= dpx) && (dy <= dpy+eps)) {
766 accepted[ico]=kFALSE;
774 // Initial value for charge ratios
775 fQrInit[0]=Float_t(fQ[fIndLocal[0][0]][0])/
776 Float_t(fQ[fIndLocal[0][0]][0]+fQ[fIndLocal[1][0]][0]);
777 fQrInit[1]=fQrInit[0];
779 if (accepted[0] && accepted[1]) {
781 fXInit[0]=0.5*(xm[0][1]+xm[0][0]);
783 fXInit[1]=0.5*(xm[0][1]+xm[1][0]);
787 chi23=CombiDoubleMathiesonFit(c);
796 } else if (accepted[0]) {
801 chi21=CombiDoubleMathiesonFit(c);
802 // Int_t ndf = fgNbins[0]+fgNbins[1]-6;
803 // Float_t prob = TMath::Prob(chi2,ndf);
804 // prob2->Fill(prob);
805 // chi2_2->Fill(chi21);
807 fprintf(stderr," chi2 %f\n",chi21);
808 if (chi21<10) Split(c);
809 } else if (accepted[1]) {
814 chi22=CombiDoubleMathiesonFit(c);
815 // Int_t ndf = fgNbins[0]+fgNbins[1]-6;
816 // Float_t prob = TMath::Prob(chi2,ndf);
817 // prob2->Fill(prob);
818 // chi2_2->Fill(chi22);
820 fprintf(stderr," chi2 %f\n",chi22);
821 if (chi22<10) Split(c);
824 if (chi21 > 10 && chi22 > 10 && chi23 > 10) {
825 // We keep only the combination found (X->cathode 2, Y->cathode 1)
826 for (Int_t ico=0; ico<2; ico++) {
828 AliMUONRawCluster cnew;
830 for (cath=0; cath<2; cath++) {
831 cnew.fX[cath]=Float_t(xm[ico][1]);
832 cnew.fY[cath]=Float_t(ym[ico][0]);
833 cnew.fZ[cath]=fZPlane;
834 cnew.fMultiplicity[cath]=c->fMultiplicity[cath];
835 for (i=0; i<fMul[cath]; i++) {
836 cnew.fIndexMap[i][cath]=c->fIndexMap[i][cath];
837 fSeg[cath]->SetPad(fIx[i][cath], fIy[i][cath]);
839 fprintf(stderr,"\nRawCluster %d cath %d\n",ico,cath);
840 fprintf(stderr,"mult_av %d\n",c->fMultiplicity[cath]);
841 FillCluster(&cnew,cath);
843 cnew.fClusterType=cnew.PhysicsContribution();
850 // +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
851 // (3') One local maximum on cathode 1 and two maxima on cathode 2
852 // +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
853 } else if (fNLocal[0]==1 && fNLocal[1]==2) {
854 Float_t xm[4][2], ym[4][2];
855 Float_t dpx, dpy, dx, dy;
856 Int_t ixm[4][2], iym[4][2];
857 Int_t isec, im1, ico;
859 // Form the 2x2 combinations
860 // 0-0, 0-1, 1-0, 1-1
862 for (im1=0; im1<2; im1++) {
863 xm[ico][0]=fX[fIndLocal[0][0]][0];
864 ym[ico][0]=fY[fIndLocal[0][0]][0];
865 xm[ico][1]=fX[fIndLocal[im1][1]][1];
866 ym[ico][1]=fY[fIndLocal[im1][1]][1];
868 ixm[ico][0]=fIx[fIndLocal[0][0]][0];
869 iym[ico][0]=fIy[fIndLocal[0][0]][0];
870 ixm[ico][1]=fIx[fIndLocal[im1][1]][1];
871 iym[ico][1]=fIy[fIndLocal[im1][1]][1];
874 // ico = 0 : first local maximum on cathodes 1 and 2
875 // ico = 1 : first local maximum on cathode 1 and second on cathode 2
877 // Analyse the combinations and keep those that are possible !
878 // For each combination check consistency in x and y
882 // In case of staggering maxima are displaced by exactly half the pad-size in y.
883 // We have to take into account the numerical precision in the consistency check;
887 for (ico=0; ico<2; ico++) {
888 accepted[ico]=kFALSE;
889 isec=fSeg[0]->Sector(ixm[ico][0], iym[ico][0]);
890 dpx=fSeg[0]->Dpx(isec)/2.;
891 dx=TMath::Abs(xm[ico][0]-xm[ico][1]);
892 isec=fSeg[1]->Sector(ixm[ico][1], iym[ico][1]);
893 dpy=fSeg[1]->Dpy(isec)/2.;
894 dy=TMath::Abs(ym[ico][0]-ym[ico][1]);
896 printf("\n %i %f %f %f %f \n", ico, ym[ico][0], ym[ico][1], dy, dpy );
897 if ((dx <= dpx) && (dy <= dpy+eps)) {
900 fprintf(stderr,"ico %d\n",ico);
904 accepted[ico]=kFALSE;
912 fQrInit[1]=Float_t(fQ[fIndLocal[0][1]][1])/
913 Float_t(fQ[fIndLocal[0][1]][1]+fQ[fIndLocal[1][1]][1]);
915 fQrInit[0]=fQrInit[1];
918 if (accepted[0] && accepted[1]) {
920 fYInit[0]=0.5*(ym[0][0]+ym[0][1]);
922 fYInit[1]=0.5*(ym[0][0]+ym[1][1]);
925 chi23=CombiDoubleMathiesonFit(c);
934 } else if (accepted[0]) {
939 chi21=CombiDoubleMathiesonFit(c);
940 // Int_t ndf = fgNbins[0]+fgNbins[1]-6;
941 // Float_t prob = TMath::Prob(chi2,ndf);
942 // prob2->Fill(prob);
943 // chi2_2->Fill(chi21);
945 fprintf(stderr," chi2 %f\n",chi21);
946 if (chi21<10) Split(c);
947 } else if (accepted[1]) {
952 chi22=CombiDoubleMathiesonFit(c);
953 // Int_t ndf = fgNbins[0]+fgNbins[1]-6;
954 // Float_t prob = TMath::Prob(chi2,ndf);
955 // prob2->Fill(prob);
956 // chi2_2->Fill(chi22);
958 fprintf(stderr," chi2 %f\n",chi22);
959 if (chi22<10) Split(c);
962 if (chi21 > 10 && chi22 > 10 && chi23 > 10) {
963 //We keep only the combination found (X->cathode 2, Y->cathode 1)
964 for (Int_t ico=0; ico<2; ico++) {
966 AliMUONRawCluster cnew;
968 for (cath=0; cath<2; cath++) {
969 cnew.fX[cath]=Float_t(xm[ico][1]);
970 cnew.fY[cath]=Float_t(ym[ico][0]);
971 cnew.fZ[cath]=fZPlane;
972 cnew.fMultiplicity[cath]=c->fMultiplicity[cath];
973 for (i=0; i<fMul[cath]; i++) {
974 cnew.fIndexMap[i][cath]=c->fIndexMap[i][cath];
975 fSeg[cath]->SetPad(fIx[i][cath], fIy[i][cath]);
977 fprintf(stderr,"\nRawCluster %d cath %d\n",ico,cath);
978 fprintf(stderr,"mult_av %d\n",c->fMultiplicity[cath]);
979 FillCluster(&cnew,cath);
981 cnew.fClusterType=cnew.PhysicsContribution();
988 // ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
989 // (4) At least three local maxima on cathode 1 or on cathode 2
990 // ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
991 } else if (fNLocal[0]>2 || fNLocal[1]>2) {
992 Int_t param = fNLocal[0]*fNLocal[1];
995 Float_t ** xm = new Float_t * [param];
996 for (ii=0; ii<param; ii++) xm[ii]=new Float_t [2];
997 Float_t ** ym = new Float_t * [param];
998 for (ii=0; ii<param; ii++) ym[ii]=new Float_t [2];
999 Int_t ** ixm = new Int_t * [param];
1000 for (ii=0; ii<param; ii++) ixm[ii]=new Int_t [2];
1001 Int_t ** iym = new Int_t * [param];
1002 for (ii=0; ii<param; ii++) iym[ii]=new Int_t [2];
1005 Float_t dpx, dpy, dx, dy;
1008 for (Int_t im1=0; im1<fNLocal[0]; im1++) {
1009 for (Int_t im2=0; im2<fNLocal[1]; im2++) {
1010 xm[ico][0]=fX[fIndLocal[im1][0]][0];
1011 ym[ico][0]=fY[fIndLocal[im1][0]][0];
1012 xm[ico][1]=fX[fIndLocal[im2][1]][1];
1013 ym[ico][1]=fY[fIndLocal[im2][1]][1];
1015 ixm[ico][0]=fIx[fIndLocal[im1][0]][0];
1016 iym[ico][0]=fIy[fIndLocal[im1][0]][0];
1017 ixm[ico][1]=fIx[fIndLocal[im2][1]][1];
1018 iym[ico][1]=fIy[fIndLocal[im2][1]][1];
1025 fprintf(stderr,"nIco %d\n",nIco);
1026 for (ico=0; ico<nIco; ico++) {
1028 fprintf(stderr,"ico = %d\n",ico);
1029 isec=fSeg[0]->Sector(ixm[ico][0], iym[ico][0]);
1030 dpx=fSeg[0]->Dpx(isec)/2.;
1031 dx=TMath::Abs(xm[ico][0]-xm[ico][1]);
1032 isec=fSeg[1]->Sector(ixm[ico][1], iym[ico][1]);
1033 dpy=fSeg[1]->Dpy(isec)/2.;
1034 dy=TMath::Abs(ym[ico][0]-ym[ico][1]);
1036 fprintf(stderr,"dx %f dpx %f dy %f dpy %f\n",dx,dpx,dy,dpy);
1037 fprintf(stderr," X %f Y %f\n",xm[ico][1],ym[ico][0]);
1039 if ((dx <= dpx) && (dy <= dpy)) {
1041 fprintf(stderr,"ok\n");
1043 AliMUONRawCluster cnew;
1044 for (cath=0; cath<2; cath++) {
1045 cnew.fX[cath]=Float_t(xm[ico][1]);
1046 cnew.fY[cath]=Float_t(ym[ico][0]);
1047 cnew.fZ[cath]=fZPlane;
1048 cnew.fMultiplicity[cath]=c->fMultiplicity[cath];
1049 for (i=0; i<fMul[cath]; i++) {
1050 cnew.fIndexMap[i][cath]=c->fIndexMap[i][cath];
1051 fSeg[cath]->SetPad(fIx[i][cath], fIy[i][cath]);
1053 FillCluster(&cnew,cath);
1055 cnew.fClusterType=cnew.PhysicsContribution();
1056 AddRawCluster(cnew);
1067 void AliMUONClusterFinderVS::FindLocalMaxima(AliMUONRawCluster* c)
1069 // Find all local maxima of a cluster
1071 printf("\n Find Local maxima !");
1075 Int_t cath, cath1; // loops over cathodes
1076 Int_t i; // loops over digits
1077 Int_t j; // loops over cathodes
1079 // Find local maxima
1081 // counters for number of local maxima
1082 fNLocal[0]=fNLocal[1]=0;
1083 // flags digits as local maximum
1084 Bool_t isLocal[100][2];
1085 for (i=0; i<100;i++) {
1086 isLocal[i][0]=isLocal[i][1]=kFALSE;
1088 // number of next neighbours and arrays to store them
1091 // loop over cathodes
1092 for (cath=0; cath<2; cath++) {
1093 // loop over cluster digits
1094 for (i=0; i<fMul[cath]; i++) {
1095 // get neighbours for that digit and assume that it is local maximum
1096 fSeg[cath]->Neighbours(fIx[i][cath], fIy[i][cath], &nn, x, y);
1097 isLocal[i][cath]=kTRUE;
1098 Int_t isec= fSeg[cath]->Sector(fIx[i][cath], fIy[i][cath]);
1099 Float_t a0 = fSeg[cath]->Dpx(isec)*fSeg[cath]->Dpy(isec);
1100 // loop over next neighbours, if at least one neighbour has higher charger assumption
1101 // digit is not local maximum
1102 for (j=0; j<nn; j++) {
1103 if (fHitMap[cath]->TestHit(x[j], y[j])==kEmpty) continue;
1104 digt=(AliMUONDigit*) fHitMap[cath]->GetHit(x[j], y[j]);
1105 isec=fSeg[cath]->Sector(x[j], y[j]);
1106 Float_t a1 = fSeg[cath]->Dpx(isec)*fSeg[cath]->Dpy(isec);
1107 if (digt->Signal()/a1 > fQ[i][cath]/a0) {
1108 isLocal[i][cath]=kFALSE;
1111 // handle special case of neighbouring pads with equal signal
1112 } else if (digt->Signal() == fQ[i][cath]) {
1113 if (fNLocal[cath]>0) {
1114 for (Int_t k=0; k<fNLocal[cath]; k++) {
1115 if (x[j]==fIx[fIndLocal[k][cath]][cath]
1116 && y[j]==fIy[fIndLocal[k][cath]][cath])
1118 isLocal[i][cath]=kFALSE;
1120 } // loop over local maxima
1121 } // are there already local maxima
1123 } // loop over next neighbours
1124 if (isLocal[i][cath]) {
1125 fIndLocal[fNLocal[cath]][cath]=i;
1128 } // loop over all digits
1129 } // loop over cathodes
1132 printf("\n Found %d %d %d %d local Maxima\n",
1133 fNLocal[0], fNLocal[1], fMul[0], fMul[1]);
1134 fprintf(stderr,"\n Cathode 1 local Maxima %d Multiplicite %d\n",fNLocal[0], fMul[0]);
1135 fprintf(stderr," Cathode 2 local Maxima %d Multiplicite %d\n",fNLocal[1], fMul[1]);
1141 if (fNLocal[1]==2 && (fNLocal[0]==1 || fNLocal[0]==0)) {
1142 Int_t iback=fNLocal[0];
1144 // Two local maxima on cathode 2 and one maximum on cathode 1
1145 // Look for local maxima considering up and down neighbours on the 1st cathode only
1147 // Loop over cluster digits
1151 for (i=0; i<fMul[cath]; i++) {
1152 isec=fSeg[cath]->Sector(fIx[i][cath],fIy[i][cath]);
1153 dpy=fSeg[cath]->Dpy(isec);
1154 dpx=fSeg[cath]->Dpx(isec);
1155 if (isLocal[i][cath]) continue;
1156 // Pad position should be consistent with position of local maxima on the opposite cathode
1157 if ((TMath::Abs(fX[i][cath]-fX[fIndLocal[0][cath1]][cath1]) > dpx/2.) &&
1158 (TMath::Abs(fX[i][cath]-fX[fIndLocal[1][cath1]][cath1]) > dpx/2.))
1161 // get neighbours for that digit and assume that it is local maximum
1162 isLocal[i][cath]=kTRUE;
1163 // compare signal to that on the two neighbours on the left and on the right
1164 // iNN counts the number of neighbours with signal, it should be 1 or 2
1168 ->FirstPad(fX[i][cath], fY[i][cath], fZPlane, 0., dpy);
1174 ix = fSeg[cath]->Ix();
1175 iy = fSeg[cath]->Iy();
1176 // skip the current pad
1177 if (iy == fIy[i][cath]) continue;
1179 if (fHitMap[cath]->TestHit(ix, iy)!=kEmpty) {
1181 digt=(AliMUONDigit*) fHitMap[cath]->GetHit(ix,iy);
1182 if (digt->Signal() > fQ[i][cath]) isLocal[i][cath]=kFALSE;
1184 } // Loop over pad neighbours in y
1185 if (isLocal[i][cath] && iNN>0) {
1186 fIndLocal[fNLocal[cath]][cath]=i;
1189 } // loop over all digits
1190 // if one additional maximum has been found we are happy
1191 // if more maxima have been found restore the previous situation
1194 "\n New search gives %d local maxima for cathode 1 \n",
1197 " %d local maxima for cathode 2 \n",
1200 if (fNLocal[cath]>2) {
1201 fNLocal[cath]=iback;
1204 } // 1,2 local maxima
1206 if (fNLocal[0]==2 && (fNLocal[1]==1 || fNLocal[1]==0)) {
1207 Int_t iback=fNLocal[1];
1209 // Two local maxima on cathode 1 and one maximum on cathode 2
1210 // Look for local maxima considering left and right neighbours on the 2nd cathode only
1213 Float_t eps = 1.e-5;
1216 // Loop over cluster digits
1217 for (i=0; i<fMul[cath]; i++) {
1218 isec=fSeg[cath]->Sector(fIx[i][cath],fIy[i][cath]);
1219 dpx=fSeg[cath]->Dpx(isec);
1220 dpy=fSeg[cath]->Dpy(isec);
1221 if (isLocal[i][cath]) continue;
1222 // Pad position should be consistent with position of local maxima on the opposite cathode
1223 if ((TMath::Abs(fY[i][cath]-fY[fIndLocal[0][cath1]][cath1]) > dpy/2.+eps) &&
1224 (TMath::Abs(fY[i][cath]-fY[fIndLocal[1][cath1]][cath1]) > dpy/2.+eps))
1228 // get neighbours for that digit and assume that it is local maximum
1229 isLocal[i][cath]=kTRUE;
1230 // compare signal to that on the two neighbours on the left and on the right
1232 // iNN counts the number of neighbours with signal, it should be 1 or 2
1235 ->FirstPad(fX[i][cath], fY[i][cath], fZPlane, dpx, 0.);
1242 ix = fSeg[cath]->Ix();
1243 iy = fSeg[cath]->Iy();
1245 // skip the current pad
1246 if (ix == fIx[i][cath]) continue;
1248 if (fHitMap[cath]->TestHit(ix, iy)!=kEmpty) {
1250 digt=(AliMUONDigit*) fHitMap[cath]->GetHit(ix,iy);
1251 if (digt->Signal() > fQ[i][cath]) isLocal[i][cath]=kFALSE;
1253 } // Loop over pad neighbours in x
1254 if (isLocal[i][cath] && iNN>0) {
1255 fIndLocal[fNLocal[cath]][cath]=i;
1258 } // loop over all digits
1259 // if one additional maximum has been found we are happy
1260 // if more maxima have been found restore the previous situation
1262 fprintf(stderr,"\n New search gives %d local maxima for cathode 1 \n",fNLocal[0]);
1263 fprintf(stderr,"\n %d local maxima for cathode 2 \n",fNLocal[1]);
1264 printf("\n New search gives %d %d \n",fNLocal[0],fNLocal[1]);
1266 if (fNLocal[cath]>2) {
1267 fNLocal[cath]=iback;
1269 } // 2,1 local maxima
1273 void AliMUONClusterFinderVS::FillCluster(AliMUONRawCluster* c, Int_t flag, Int_t cath)
1276 // Completes cluster information starting from list of digits
1283 c->fPeakSignal[cath]=c->fPeakSignal[0];
1285 c->fPeakSignal[cath]=0;
1296 fprintf(stderr,"\n fPeakSignal %d\n",c->fPeakSignal[cath]);
1297 for (Int_t i=0; i<c->fMultiplicity[cath]; i++)
1299 dig= fInput->Digit(cath,c->fIndexMap[i][cath]);
1300 ix=dig->PadX()+c->fOffsetMap[i][cath];
1302 Int_t q=dig->Signal();
1303 if (!flag) q=Int_t(q*c->fContMap[i][cath]);
1304 // fprintf(stderr,"q %d c->fPeakSignal[ %d ] %d\n",q,cath,c->fPeakSignal[cath]);
1305 if (dig->Physics() >= dig->Signal()) {
1306 c->fPhysicsMap[i]=2;
1307 } else if (dig->Physics() == 0) {
1308 c->fPhysicsMap[i]=0;
1309 } else c->fPhysicsMap[i]=1;
1313 fprintf(stderr,"q %d c->fPeakSignal[cath] %d\n",q,c->fPeakSignal[cath]);
1314 // peak signal and track list
1315 if (q>c->fPeakSignal[cath]) {
1316 c->fPeakSignal[cath]=q;
1317 c->fTracks[0]=dig->Hit();
1318 c->fTracks[1]=dig->Track(0);
1319 c->fTracks[2]=dig->Track(1);
1320 // fprintf(stderr," c->fTracks[0] %d c->fTracks[1] %d\n",dig->fHit,dig->fTracks[0]);
1324 fSeg[cath]->GetPadC(ix, iy, x, y, z);
1329 } // loop over digits
1331 fprintf(stderr," fin du cluster c\n");
1335 c->fX[cath]/=c->fQ[cath];
1337 c->fX[cath]=fSeg[cath]->GetAnod(c->fX[cath]);
1338 c->fY[cath]/=c->fQ[cath];
1340 // apply correction to the coordinate along the anode wire
1344 fSeg[cath]->GetPadI(x, y, fZPlane, ix, iy);
1345 fSeg[cath]->GetPadC(ix, iy, x, y, z);
1346 Int_t isec=fSeg[cath]->Sector(ix,iy);
1347 TF1* cogCorr = fSeg[cath]->CorrFunc(isec-1);
1350 Float_t yOnPad=(c->fY[cath]-y)/fSeg[cath]->Dpy(isec);
1351 c->fY[cath]=c->fY[cath]-cogCorr->Eval(yOnPad, 0, 0);
1356 void AliMUONClusterFinderVS::FillCluster(AliMUONRawCluster* c, Int_t cath)
1359 // Completes cluster information starting from list of digits
1369 Float_t xpad, ypad, zpad;
1372 for (Int_t i=0; i<c->fMultiplicity[cath]; i++)
1374 dig = fInput->Digit(cath,c->fIndexMap[i][cath]);
1376 GetPadC(dig->PadX(),dig->PadY(),xpad,ypad, zpad);
1378 fprintf(stderr,"x %f y %f cx %f cy %f\n",xpad,ypad,c->fX[0],c->fY[0]);
1379 dx = xpad - c->fX[0];
1380 dy = ypad - c->fY[0];
1381 dr = TMath::Sqrt(dx*dx+dy*dy);
1386 fprintf(stderr," dr %f\n",dr);
1387 Int_t q=dig->Signal();
1388 if (dig->Physics() >= dig->Signal()) {
1389 c->fPhysicsMap[i]=2;
1390 } else if (dig->Physics() == 0) {
1391 c->fPhysicsMap[i]=0;
1392 } else c->fPhysicsMap[i]=1;
1393 c->fPeakSignal[cath]=q;
1394 c->fTracks[0]=dig->Hit();
1395 c->fTracks[1]=dig->Track(0);
1396 c->fTracks[2]=dig->Track(1);
1398 fprintf(stderr," c->fTracks[0] %d c->fTracks[1] %d\n",dig->Hit(),
1402 } // loop over digits
1404 // apply correction to the coordinate along the anode wire
1406 c->fX[cath]=fSeg[cath]->GetAnod(c->fX[cath]);
1409 void AliMUONClusterFinderVS::FindCluster(Int_t i, Int_t j, Int_t cath, AliMUONRawCluster &c){
1413 // Find a super cluster on both cathodes
1416 // Add i,j as element of the cluster
1419 Int_t idx = fHitMap[cath]->GetHitIndex(i,j);
1420 AliMUONDigit* dig = (AliMUONDigit*) fHitMap[cath]->GetHit(i,j);
1421 Int_t q=dig->Signal();
1422 Int_t theX=dig->PadX();
1423 Int_t theY=dig->PadY();
1425 if (q > TMath::Abs(c.fPeakSignal[0]) && q > TMath::Abs(c.fPeakSignal[1])) {
1426 c.fPeakSignal[cath]=q;
1427 c.fTracks[0]=dig->Hit();
1428 c.fTracks[1]=dig->Track(0);
1429 c.fTracks[2]=dig->Track(1);
1433 // Make sure that list of digits is ordered
1435 Int_t mu=c.fMultiplicity[cath];
1436 c.fIndexMap[mu][cath]=idx;
1438 if (dig->Physics() >= dig->Signal()) {
1439 c.fPhysicsMap[mu]=2;
1440 } else if (dig->Physics() == 0) {
1441 c.fPhysicsMap[mu]=0;
1442 } else c.fPhysicsMap[mu]=1;
1446 for (Int_t ind = mu-1; ind >= 0; ind--) {
1447 Int_t ist=(c.fIndexMap)[ind][cath];
1448 Int_t ql=fInput->Digit(cath, ist)->Signal();
1449 Int_t ix=fInput->Digit(cath, ist)->PadX();
1450 Int_t iy=fInput->Digit(cath, ist)->PadY();
1452 if (q>ql || (q==ql && theX > ix && theY < iy)) {
1453 c.fIndexMap[ind][cath]=idx;
1454 c.fIndexMap[ind+1][cath]=ist;
1462 c.fMultiplicity[cath]++;
1463 if (c.fMultiplicity[cath] >= 50 ) {
1464 printf("FindCluster - multiplicity >50 %d \n",c.fMultiplicity[0]);
1465 c.fMultiplicity[cath]=49;
1468 // Prepare center of gravity calculation
1470 fSeg[cath]->GetPadC(i, j, x, y, z);
1476 // Flag hit as "taken"
1477 fHitMap[cath]->FlagHit(i,j);
1479 // Now look recursively for all neighbours and pad hit on opposite cathode
1481 // Loop over neighbours
1485 Int_t xList[10], yList[10];
1486 fSeg[cath]->Neighbours(i,j,&nn,xList,yList);
1487 for (Int_t in=0; in<nn; in++) {
1491 if (fHitMap[cath]->TestHit(ix,iy)==kUnused) {
1493 printf("\n Neighbours %d %d %d", cath, ix, iy);
1494 FindCluster(ix, iy, cath, c);
1499 Int_t iXopp[50], iYopp[50];
1501 // Neighbours on opposite cathode
1502 // Take into account that several pads can overlap with the present pad
1503 Int_t isec=fSeg[cath]->Sector(i,j);
1509 dx = (fSeg[cath]->Dpx(isec))/2.;
1514 dy = (fSeg[cath]->Dpy(isec))/2;
1516 // loop over pad neighbours on opposite cathode
1517 for (fSeg[iop]->FirstPad(x, y, fZPlane, dx, dy);
1518 fSeg[iop]->MorePads();
1519 fSeg[iop]->NextPad())
1522 ix = fSeg[iop]->Ix(); iy = fSeg[iop]->Iy();
1523 if (fDebugLevel > 1)
1524 printf("\n ix, iy: %f %f %f %d %d %d", x,y,z,ix, iy, fSector);
1525 if (fHitMap[iop]->TestHit(ix,iy)==kUnused){
1528 if (fDebugLevel > 1)
1529 printf("\n Opposite %d %d %d", iop, ix, iy);
1532 } // Loop over pad neighbours
1533 // This had to go outside the loop since recursive calls inside the iterator are not possible
1536 for (jopp=0; jopp<nOpp; jopp++) {
1537 if (fHitMap[iop]->TestHit(iXopp[jopp],iYopp[jopp]) == kUnused)
1538 FindCluster(iXopp[jopp], iYopp[jopp], iop, c);
1542 //_____________________________________________________________________________
1544 void AliMUONClusterFinderVS::FindRawClusters()
1547 // MUON cluster finder from digits -- finds neighbours on both cathodes and
1548 // fills the tree with raw clusters
1551 // Return if no input datad available
1552 if (!fInput->NDigits(0) && !fInput->NDigits(1)) return;
1554 fSeg[0] = fInput->Segmentation(0);
1555 fSeg[1] = fInput->Segmentation(1);
1557 fHitMap[0] = new AliMUONHitMapA1(fSeg[0], fInput->Digits(0));
1558 fHitMap[1] = new AliMUONHitMapA1(fSeg[1], fInput->Digits(1));
1566 fHitMap[0]->FillHits();
1567 fHitMap[1]->FillHits();
1569 // Outer Loop over Cathodes
1570 for (cath=0; cath<2; cath++) {
1571 for (ndig=0; ndig<fInput->NDigits(cath); ndig++) {
1572 dig = fInput->Digit(cath, ndig);
1573 Int_t i=dig->PadX();
1574 Int_t j=dig->PadY();
1575 if (fHitMap[cath]->TestHit(i,j)==kUsed ||fHitMap[0]->TestHit(i,j)==kEmpty) {
1580 fprintf(stderr,"\n CATHODE %d CLUSTER %d\n",cath,ncls);
1581 AliMUONRawCluster c;
1582 c.fMultiplicity[0]=0;
1583 c.fMultiplicity[1]=0;
1584 c.fPeakSignal[cath]=dig->Signal();
1585 c.fTracks[0]=dig->Hit();
1586 c.fTracks[1]=dig->Track(0);
1587 c.fTracks[2]=dig->Track(1);
1588 // tag the beginning of cluster list in a raw cluster
1591 fSeg[cath]->GetPadC(i,j,xcu, ycu, fZPlane);
1592 fSector= fSeg[cath]->Sector(i,j)/100;
1594 printf("\n New Seed %d %d ", i,j);
1596 FindCluster(i,j,cath,c);
1597 // ^^^^^^^^^^^^^^^^^^^^^^^^
1598 // center of gravity
1601 c.fX[0]=fSeg[0]->GetAnod(c.fX[0]);
1605 c.fX[1]=fSeg[0]->GetAnod(c.fX[1]);
1612 fprintf(stderr,"\n Cathode 1 multiplicite %d X(CG) %f Y(CG) %f\n",
1613 c.fMultiplicity[0],c.fX[0],c.fY[0]);
1614 fprintf(stderr," Cathode 2 multiplicite %d X(CG) %f Y(CG) %f\n",
1615 c.fMultiplicity[1],c.fX[1],c.fY[1]);
1617 // Analyse cluster and decluster if necessary
1620 c.fNcluster[1]=fNRawClusters;
1621 c.fClusterType=c.PhysicsContribution();
1628 // reset Cluster object
1629 { // begin local scope
1630 for (int k=0;k<c.fMultiplicity[0];k++) c.fIndexMap[k][0]=0;
1631 } // end local scope
1633 { // begin local scope
1634 for (int k=0;k<c.fMultiplicity[1];k++) c.fIndexMap[k][1]=0;
1635 } // end local scope
1637 c.fMultiplicity[0]=c.fMultiplicity[0]=0;
1641 } // end loop cathodes
1646 Float_t AliMUONClusterFinderVS::SingleMathiesonFit(AliMUONRawCluster *c, Int_t cath)
1648 // Performs a single Mathieson fit on one cathode
1650 AliMUONClusterInput& clusterInput = *(AliMUONClusterInput::Instance());
1652 clusterInput.Fitter()->SetFCN(fcnS1);
1653 clusterInput.Fitter()->mninit(2,10,7);
1654 Double_t arglist[20];
1657 // Set starting values
1658 static Double_t vstart[2];
1663 // lower and upper limits
1664 static Double_t lower[2], upper[2];
1666 fSeg[cath]->GetPadI(c->fX[cath], c->fY[cath], fZPlane, ix, iy);
1667 Int_t isec=fSeg[cath]->Sector(ix, iy);
1668 lower[0]=vstart[0]-fSeg[cath]->Dpx(isec)/2;
1669 lower[1]=vstart[1]-fSeg[cath]->Dpy(isec)/2;
1671 upper[0]=lower[0]+fSeg[cath]->Dpx(isec);
1672 upper[1]=lower[1]+fSeg[cath]->Dpy(isec);
1675 static Double_t step[2]={0.0005, 0.0005};
1677 clusterInput.Fitter()->mnparm(0,"x1",vstart[0],step[0],lower[0],upper[0],ierflag);
1678 clusterInput.Fitter()->mnparm(1,"y1",vstart[1],step[1],lower[1],upper[1],ierflag);
1679 // ready for minimisation
1680 clusterInput.Fitter()->SetPrintLevel(-1+fDebugLevel);
1682 clusterInput.Fitter()->mnexcm("SET NOW", arglist, 0, ierflag);
1683 clusterInput.Fitter()->mnexcm("SET OUT", arglist, 0, ierflag);
1687 clusterInput.Fitter()->mnexcm("SET NOGR", arglist, 0, ierflag);
1688 clusterInput.Fitter()->mnexcm("MIGRAD", arglist, 0, ierflag);
1689 clusterInput.Fitter()->mnexcm("EXIT" , arglist, 0, ierflag);
1690 Double_t fmin, fedm, errdef;
1691 Int_t npari, nparx, istat;
1693 clusterInput.Fitter()->mnstat(fmin, fedm, errdef, npari, nparx, istat);
1697 // Get fitted parameters
1698 Double_t xrec, yrec;
1700 Double_t epxz, b1, b2;
1702 clusterInput.Fitter()->mnpout(0, chname, xrec, epxz, b1, b2, ierflg);
1703 clusterInput.Fitter()->mnpout(1, chname, yrec, epxz, b1, b2, ierflg);
1709 Float_t AliMUONClusterFinderVS::CombiSingleMathiesonFit(AliMUONRawCluster *c)
1711 // Perform combined Mathieson fit on both cathode planes
1713 AliMUONClusterInput& clusterInput = *(AliMUONClusterInput::Instance());
1714 clusterInput.Fitter()->SetFCN(fcnCombiS1);
1715 clusterInput.Fitter()->mninit(2,10,7);
1716 Double_t arglist[20];
1719 static Double_t vstart[2];
1720 vstart[0]=fXInit[0];
1721 vstart[1]=fYInit[0];
1724 // lower and upper limits
1725 static Float_t lower[2], upper[2];
1727 fSeg[0]->GetPadI(fXInit[0], fYInit[0], fZPlane, ix, iy);
1728 isec=fSeg[0]->Sector(ix, iy);
1729 Float_t dpy=fSeg[0]->Dpy(isec);
1730 fSeg[1]->GetPadI(fXInit[0], fYInit[0], fZPlane, ix, iy);
1731 isec=fSeg[1]->Sector(ix, iy);
1732 Float_t dpx=fSeg[1]->Dpx(isec);
1735 Float_t xdum, ydum, zdum;
1737 // Find save upper and lower limits
1741 for (fSeg[1]->FirstPad(fXInit[0], fYInit[0], fZPlane, dpx, 0.);
1742 fSeg[1]->MorePads(); fSeg[1]->NextPad())
1744 ix=fSeg[1]->Ix(); iy=fSeg[1]->Iy();
1745 fSeg[1]->GetPadC(ix,iy, upper[0], ydum, zdum);
1746 if (icount ==0) lower[0]=upper[0];
1750 if (lower[0]>upper[0]) {xdum=lower[0]; lower[0]=upper[0]; upper[0]=xdum;}
1754 printf("\n single y %f %f", fXInit[0], fYInit[0]);
1756 for (fSeg[0]->FirstPad(fXInit[0], fYInit[0], fZPlane, 0., dpy);
1757 fSeg[0]->MorePads(); fSeg[0]->NextPad())
1759 ix=fSeg[0]->Ix(); iy=fSeg[0]->Iy();
1760 fSeg[0]->GetPadC(ix,iy,xdum,upper[1],zdum);
1761 if (icount ==0) lower[1]=upper[1];
1764 printf("\n upper lower %d %f %f", icount, upper[1], lower[1]);
1767 if (lower[1]>upper[1]) {xdum=lower[1]; lower[1]=upper[1]; upper[1]=xdum;}
1770 static Double_t step[2]={0.00001, 0.0001};
1772 clusterInput.Fitter()->mnparm(0,"x1",vstart[0],step[0],lower[0],upper[0],ierflag);
1773 clusterInput.Fitter()->mnparm(1,"y1",vstart[1],step[1],lower[1],upper[1],ierflag);
1774 // ready for minimisation
1775 clusterInput.Fitter()->SetPrintLevel(-1+fDebugLevel);
1777 clusterInput.Fitter()->mnexcm("SET NOW", arglist, 0, ierflag);
1778 clusterInput.Fitter()->mnexcm("SET OUT", arglist, 0, ierflag);
1782 clusterInput.Fitter()->mnexcm("SET NOGR", arglist, 0, ierflag);
1783 clusterInput.Fitter()->mnexcm("MIGRAD", arglist, 0, ierflag);
1784 clusterInput.Fitter()->mnexcm("EXIT" , arglist, 0, ierflag);
1785 Double_t fmin, fedm, errdef;
1786 Int_t npari, nparx, istat;
1788 clusterInput.Fitter()->mnstat(fmin, fedm, errdef, npari, nparx, istat);
1792 // Get fitted parameters
1793 Double_t xrec, yrec;
1795 Double_t epxz, b1, b2;
1797 clusterInput.Fitter()->mnpout(0, chname, xrec, epxz, b1, b2, ierflg);
1798 clusterInput.Fitter()->mnpout(1, chname, yrec, epxz, b1, b2, ierflg);
1804 Bool_t AliMUONClusterFinderVS::DoubleMathiesonFit(AliMUONRawCluster *c, Int_t cath)
1806 // Performs a double Mathieson fit on one cathode
1810 // Initialise global variables for fit
1811 AliMUONClusterInput& clusterInput = *(AliMUONClusterInput::Instance());
1812 clusterInput.Fitter()->SetFCN(fcnS2);
1813 clusterInput.Fitter()->mninit(5,10,7);
1814 Double_t arglist[20];
1817 // Set starting values
1818 static Double_t vstart[5];
1819 vstart[0]=fX[fIndLocal[0][cath]][cath];
1820 vstart[1]=fY[fIndLocal[0][cath]][cath];
1821 vstart[2]=fX[fIndLocal[1][cath]][cath];
1822 vstart[3]=fY[fIndLocal[1][cath]][cath];
1823 vstart[4]=Float_t(fQ[fIndLocal[0][cath]][cath])/
1824 Float_t(fQ[fIndLocal[0][cath]][cath]+fQ[fIndLocal[1][cath]][cath]);
1825 // lower and upper limits
1826 static Float_t lower[5], upper[5];
1827 Int_t isec=fSeg[cath]->Sector(fIx[fIndLocal[0][cath]][cath], fIy[fIndLocal[0][cath]][cath]);
1828 lower[0]=vstart[0]-fSeg[cath]->Dpx(isec);
1829 lower[1]=vstart[1]-fSeg[cath]->Dpy(isec);
1831 upper[0]=lower[0]+2.*fSeg[cath]->Dpx(isec);
1832 upper[1]=lower[1]+2.*fSeg[cath]->Dpy(isec);
1834 isec=fSeg[cath]->Sector(fIx[fIndLocal[1][cath]][cath], fIy[fIndLocal[1][cath]][cath]);
1835 lower[2]=vstart[2]-fSeg[cath]->Dpx(isec)/2;
1836 lower[3]=vstart[3]-fSeg[cath]->Dpy(isec)/2;
1838 upper[2]=lower[2]+fSeg[cath]->Dpx(isec);
1839 upper[3]=lower[3]+fSeg[cath]->Dpy(isec);
1844 static Double_t step[5]={0.0005, 0.0005, 0.0005, 0.0005, 0.0001};
1846 clusterInput.Fitter()->mnparm(0,"x1",vstart[0],step[0],lower[0],upper[0],ierflag);
1847 clusterInput.Fitter()->mnparm(1,"y1",vstart[1],step[1],lower[1],upper[1],ierflag);
1848 clusterInput.Fitter()->mnparm(2,"x2",vstart[2],step[2],lower[2],upper[2],ierflag);
1849 clusterInput.Fitter()->mnparm(3,"y2",vstart[3],step[3],lower[3],upper[3],ierflag);
1850 clusterInput.Fitter()->mnparm(4,"a0",vstart[4],step[4],lower[4],upper[4],ierflag);
1851 // ready for minimisation
1852 clusterInput.Fitter()->SetPrintLevel(-1+fDebugLevel);
1854 clusterInput.Fitter()->mnexcm("SET NOW", arglist, 0, ierflag);
1855 clusterInput.Fitter()->mnexcm("SET OUT", arglist, 0, ierflag);
1859 clusterInput.Fitter()->mnexcm("SET NOGR", arglist, 0, ierflag);
1860 clusterInput.Fitter()->mnexcm("MIGRAD", arglist, 0, ierflag);
1861 clusterInput.Fitter()->mnexcm("EXIT" , arglist, 0, ierflag);
1862 // Get fitted parameters
1863 Double_t xrec[2], yrec[2], qfrac;
1865 Double_t epxz, b1, b2;
1867 clusterInput.Fitter()->mnpout(0, chname, xrec[0], epxz, b1, b2, ierflg);
1868 clusterInput.Fitter()->mnpout(1, chname, yrec[0], epxz, b1, b2, ierflg);
1869 clusterInput.Fitter()->mnpout(2, chname, xrec[1], epxz, b1, b2, ierflg);
1870 clusterInput.Fitter()->mnpout(3, chname, yrec[1], epxz, b1, b2, ierflg);
1871 clusterInput.Fitter()->mnpout(4, chname, qfrac, epxz, b1, b2, ierflg);
1873 Double_t fmin, fedm, errdef;
1874 Int_t npari, nparx, istat;
1876 clusterInput.Fitter()->mnstat(fmin, fedm, errdef, npari, nparx, istat);
1881 Float_t AliMUONClusterFinderVS::CombiDoubleMathiesonFit(AliMUONRawCluster *c)
1884 // Perform combined double Mathieson fit on both cathode planes
1886 AliMUONClusterInput& clusterInput = *(AliMUONClusterInput::Instance());
1887 clusterInput.Fitter()->SetFCN(fcnCombiS2);
1888 clusterInput.Fitter()->mninit(6,10,7);
1889 Double_t arglist[20];
1892 // Set starting values
1893 static Double_t vstart[6];
1894 vstart[0]=fXInit[0];
1895 vstart[1]=fYInit[0];
1896 vstart[2]=fXInit[1];
1897 vstart[3]=fYInit[1];
1898 vstart[4]=fQrInit[0];
1899 vstart[5]=fQrInit[1];
1900 // lower and upper limits
1901 static Float_t lower[6], upper[6];
1905 fSeg[1]->GetPadI(fXInit[0], fYInit[0], fZPlane, ix, iy);
1906 isec=fSeg[1]->Sector(ix, iy);
1907 dpx=fSeg[1]->Dpx(isec);
1909 fSeg[0]->GetPadI(fXInit[0], fYInit[0], fZPlane, ix, iy);
1910 isec=fSeg[0]->Sector(ix, iy);
1911 dpy=fSeg[0]->Dpy(isec);
1915 Float_t xdum, ydum, zdum;
1917 printf("\n Cluster Finder: %f %f %f %f ", fXInit[0], fXInit[1],fYInit[0], fYInit[1] );
1919 // Find save upper and lower limits
1922 for (fSeg[1]->FirstPad(fXInit[0], fYInit[0], fZPlane, dpx, 0.);
1923 fSeg[1]->MorePads(); fSeg[1]->NextPad())
1925 ix=fSeg[1]->Ix(); iy=fSeg[1]->Iy();
1926 // if (fHitMap[1]->TestHit(ix, iy) == kEmpty) continue;
1927 fSeg[1]->GetPadC(ix,iy,upper[0],ydum,zdum);
1928 if (icount ==0) lower[0]=upper[0];
1931 if (lower[0]>upper[0]) {xdum=lower[0]; lower[0]=upper[0]; upper[0]=xdum;}
1932 // vstart[0] = 0.5*(lower[0]+upper[0]);
1937 for (fSeg[0]->FirstPad(fXInit[0], fYInit[0], fZPlane, 0., dpy);
1938 fSeg[0]->MorePads(); fSeg[0]->NextPad())
1940 ix=fSeg[0]->Ix(); iy=fSeg[0]->Iy();
1941 // if (fHitMap[0]->TestHit(ix, iy) == kEmpty) continue;
1942 fSeg[0]->GetPadC(ix,iy,xdum,upper[1],zdum);
1943 if (icount ==0) lower[1]=upper[1];
1947 if (lower[1]>upper[1]) {xdum=lower[1]; lower[1]=upper[1]; upper[1]=xdum;}
1948 // vstart[1] = 0.5*(lower[1]+upper[1]);
1951 fSeg[1]->GetPadI(fXInit[1], fYInit[1], fZPlane, ix, iy);
1952 isec=fSeg[1]->Sector(ix, iy);
1953 dpx=fSeg[1]->Dpx(isec);
1954 fSeg[0]->GetPadI(fXInit[1], fYInit[1], fZPlane, ix, iy);
1955 isec=fSeg[0]->Sector(ix, iy);
1956 dpy=fSeg[0]->Dpy(isec);
1959 // Find save upper and lower limits
1963 for (fSeg[1]->FirstPad(fXInit[1], fYInit[1], fZPlane, dpx, 0);
1964 fSeg[1]->MorePads(); fSeg[1]->NextPad())
1966 ix=fSeg[1]->Ix(); iy=fSeg[1]->Iy();
1967 // if (fHitMap[1]->TestHit(ix, iy) == kEmpty) continue;
1968 fSeg[1]->GetPadC(ix,iy,upper[2],ydum,zdum);
1969 if (icount ==0) lower[2]=upper[2];
1972 if (lower[2]>upper[2]) {xdum=lower[2]; lower[2]=upper[2]; upper[2]=xdum;}
1973 // vstart[2] = 0.5*(lower[2]+upper[2]);
1977 for (fSeg[0]->FirstPad(fXInit[1], fYInit[1], fZPlane, 0, dpy);
1978 fSeg[0]-> MorePads(); fSeg[0]->NextPad())
1980 ix=fSeg[0]->Ix(); iy=fSeg[0]->Iy();
1981 // if (fHitMap[0]->TestHit(ix, iy) != kEmpty) continue;
1983 fSeg[0]->GetPadC(ix,iy,xdum,upper[3],zdum);
1984 if (icount ==0) lower[3]=upper[3];
1988 if (lower[3]>upper[3]) {xdum=lower[3]; lower[3]=upper[3]; upper[3]=xdum;}
1990 // vstart[3] = 0.5*(lower[3]+upper[3]);
1998 static Double_t step[6]={0.0005, 0.0005, 0.0005, 0.0005, 0.001, 0.001};
1999 clusterInput.Fitter()->mnparm(0,"x1",vstart[0],step[0],lower[0],upper[0],ierflag);
2000 clusterInput.Fitter()->mnparm(1,"y1",vstart[1],step[1],lower[1],upper[1],ierflag);
2001 clusterInput.Fitter()->mnparm(2,"x2",vstart[2],step[2],lower[2],upper[2],ierflag);
2002 clusterInput.Fitter()->mnparm(3,"y2",vstart[3],step[3],lower[3],upper[3],ierflag);
2003 clusterInput.Fitter()->mnparm(4,"a0",vstart[4],step[4],lower[4],upper[4],ierflag);
2004 clusterInput.Fitter()->mnparm(5,"a1",vstart[5],step[5],lower[5],upper[5],ierflag);
2005 // ready for minimisation
2006 clusterInput.Fitter()->SetPrintLevel(-1+fDebugLevel);
2008 clusterInput.Fitter()->mnexcm("SET NOW", arglist, 0, ierflag);
2009 clusterInput.Fitter()->mnexcm("SET OUT", arglist, 0, ierflag);
2013 clusterInput.Fitter()->mnexcm("SET NOGR", arglist, 0, ierflag);
2014 clusterInput.Fitter()->mnexcm("MIGRAD", arglist, 0, ierflag);
2015 clusterInput.Fitter()->mnexcm("EXIT" , arglist, 0, ierflag);
2016 // Get fitted parameters
2018 Double_t epxz, b1, b2;
2020 clusterInput.Fitter()->mnpout(0, chname, fXFit[0], epxz, b1, b2, ierflg);
2021 clusterInput.Fitter()->mnpout(1, chname, fYFit[0], epxz, b1, b2, ierflg);
2022 clusterInput.Fitter()->mnpout(2, chname, fXFit[1], epxz, b1, b2, ierflg);
2023 clusterInput.Fitter()->mnpout(3, chname, fYFit[1], epxz, b1, b2, ierflg);
2024 clusterInput.Fitter()->mnpout(4, chname, fQrFit[0], epxz, b1, b2, ierflg);
2025 clusterInput.Fitter()->mnpout(5, chname, fQrFit[1], epxz, b1, b2, ierflg);
2027 Double_t fmin, fedm, errdef;
2028 Int_t npari, nparx, istat;
2030 clusterInput.Fitter()->mnstat(fmin, fedm, errdef, npari, nparx, istat);
2038 void AliMUONClusterFinderVS::Split(AliMUONRawCluster* c)
2041 // One cluster for each maximum
2044 AliMUONClusterInput& clusterInput = *(AliMUONClusterInput::Instance());
2045 for (j=0; j<2; j++) {
2046 AliMUONRawCluster cnew;
2047 cnew.fGhost=c->fGhost;
2048 for (cath=0; cath<2; cath++) {
2049 cnew.fChi2[cath]=fChi2[0];
2050 // ?? why not cnew.fChi2[cath]=fChi2[cath];
2053 cnew.fNcluster[0]=-1;
2054 cnew.fNcluster[1]=fNRawClusters;
2056 cnew.fNcluster[0]=fNPeaks;
2057 cnew.fNcluster[1]=0;
2059 cnew.fMultiplicity[cath]=0;
2060 cnew.fX[cath]=Float_t(fXFit[j]);
2061 cnew.fY[cath]=Float_t(fYFit[j]);
2062 cnew.fZ[cath]=fZPlane;
2064 cnew.fQ[cath]=Int_t(clusterInput.TotalCharge(cath)*fQrFit[cath]);
2066 cnew.fQ[cath]=Int_t(clusterInput.TotalCharge(cath)*(1-fQrFit[cath]));
2068 fSeg[cath]->SetHit(fXFit[j],fYFit[j],fZPlane);
2069 for (i=0; i<fMul[cath]; i++) {
2070 cnew.fIndexMap[cnew.fMultiplicity[cath]][cath]=
2071 c->fIndexMap[i][cath];
2072 fSeg[cath]->SetPad(fIx[i][cath], fIy[i][cath]);
2073 Float_t q1=fInput->Response()->IntXY(fSeg[cath]);
2074 cnew.fContMap[i][cath]
2075 =(q1*Float_t(cnew.fQ[cath]))/Float_t(fQ[i][cath]);
2076 cnew.fMultiplicity[cath]++;
2078 FillCluster(&cnew,0,cath);
2081 cnew.fClusterType=cnew.PhysicsContribution();
2082 if (cnew.fQ[0]>0 && cnew.fQ[1]>0) AddRawCluster(cnew);
2089 // Minimisation functions
2091 void fcnS1(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag)
2093 AliMUONClusterInput& clusterInput = *(AliMUONClusterInput::Instance());
2100 for (i=0; i<clusterInput.Nmul(0); i++) {
2101 Float_t q0=clusterInput.Charge(i,0);
2102 Float_t q1=clusterInput.DiscrChargeS1(i,par);
2111 void fcnCombiS1(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag)
2113 AliMUONClusterInput& clusterInput = *(AliMUONClusterInput::Instance());
2120 for (cath=0; cath<2; cath++) {
2121 for (i=0; i<clusterInput.Nmul(cath); i++) {
2122 Float_t q0=clusterInput.Charge(i,cath);
2123 Float_t q1=clusterInput.DiscrChargeCombiS1(i,par,cath);
2134 void fcnS2(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag)
2136 AliMUONClusterInput& clusterInput = *(AliMUONClusterInput::Instance());
2143 for (i=0; i<clusterInput.Nmul(0); i++) {
2145 Float_t q0=clusterInput.Charge(i,0);
2146 Float_t q1=clusterInput.DiscrChargeS2(i,par);
2156 void fcnCombiS2(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag)
2158 AliMUONClusterInput& clusterInput = *(AliMUONClusterInput::Instance());
2164 for (cath=0; cath<2; cath++) {
2165 for (i=0; i<clusterInput.Nmul(cath); i++) {
2166 Float_t q0=clusterInput.Charge(i,cath);
2167 Float_t q1=clusterInput.DiscrChargeCombiS2(i,par,cath);
2177 void AliMUONClusterFinderVS::AddRawCluster(const AliMUONRawCluster c)
2180 // Add a raw cluster copy to the list
2182 AliMUON *pMUON=(AliMUON*)gAlice->GetModule("MUON");
2183 pMUON->AddRawCluster(fInput->Chamber(),c);
2186 fprintf(stderr,"\nfNRawClusters %d\n",fNRawClusters);
2189 Bool_t AliMUONClusterFinderVS::TestTrack(Int_t t) {
2190 // Test if track was user selected
2191 if (fTrack[0]==-1 || fTrack[1]==-1) {
2193 } else if (t==fTrack[0] || t==fTrack[1]) {
2200 AliMUONClusterFinderVS& AliMUONClusterFinderVS
2201 ::operator = (const AliMUONClusterFinderVS& rhs)
2203 // Dummy assignment operator