]> git.uio.no Git - u/mrichter/AliRoot.git/blame_incremental - MUON/AliMUONClusterFinderVS.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame_incremental
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (non-incremental) | history | HEAD
Bug in decision on splitting corrected.
[u/mrichter/AliRoot.git] / MUON / AliMUONClusterFinderVS.cxx
... / ...
CommitLineData
1/**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
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 **************************************************************************/
15/*
16$Log$
17Revision 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
21
22Revision 1.19 2001/03/20 13:32:10 egangler
23Code introduced to remove ghosts with the charge correlation between the 2
24cathods. A chi2 is performed for the 2 possibilities.
25If one gets good chi2 (with respect to the fGhostChi2Cut parameter) and the
26other wrong chi2, the ambiguity is solved
27If both gets good or both bad chi2, no choice is made
28By default the fGhostChi2Cut parameter is set to solve about 70% of ghost
29problems with about 2% errors, with the current version of the code.
30
31Implementation :
32fGhostChi2Cut is in AliMUONClusterFinderVS, with setters and getters.
33a fDebugLevel was also introduced to switch off some of the output.
34When an ambiguity is detected and not solved, the fGhost word in
35AliMUONRawCluster is set to 1 or 2, depending whether both charge chi2 are
36good or bad.
37a DumpIndex method was also added in AliMUONRawCluster to produce a printout
38of digit indexes.
39
40User incidences :
41By default, the code makes ghost check. If you want previous behaviour,
42put in MUONrawclusters the value of SetGhostChi2Cut to infinity (1e.6) is
43sufficient.
44
45Revision 1.18 2001/01/26 21:37:53 morsch
46Use access functions to AliMUONDigit member data.
47
48Revision 1.17 2001/01/23 18:58:19 hristov
49Initialisation of some pointers
50
51Revision 1.16 2000/12/21 23:27:30 morsch
52Error in argument list of AddRawCluster corrected.
53
54Revision 1.15 2000/12/21 22:14:38 morsch
55Clean-up of coding rule violations.
56
57Revision 1.14 2000/10/23 16:03:45 morsch
58Correct z-position of all clusters created "on the flight".
59
60Revision 1.13 2000/10/23 13:38:23 morsch
61Set correct z-coordinate when cluster is split.
62
63Revision 1.12 2000/10/18 11:42:06 morsch
64- AliMUONRawCluster contains z-position.
65- Some clean-up of useless print statements during initialisations.
66
67Revision 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.
71
72Revision 1.10 2000/10/03 13:51:57 egangler
73Removal of useless dependencies via forward declarations
74
75Revision 1.9 2000/10/02 16:58:29 egangler
76Cleaning of the code :
77-> coding conventions
78-> void Streamers
79-> some useless includes removed or replaced by "class" statement
80
81Revision 1.8 2000/07/03 11:54:57 morsch
82AliMUONSegmentation and AliMUONHitMap have been replaced by AliSegmentation and AliHitMap in STEER
83The methods GetPadIxy and GetPadXxy of AliMUONSegmentation have changed name to GetPadI and GetPadC.
84
85Revision 1.7 2000/06/28 15:16:35 morsch
86(1) Client code adapted to new method signatures in AliMUONSegmentation (see comments there)
87to allow development of slat-muon chamber simulation and reconstruction code in the MUON
88framework. The changes should have no side effects (mostly dummy arguments).
89(2) Hit disintegration uses 3-dim hit coordinates to allow simulation
90of chambers with overlapping modules (MakePadHits, Disintegration).
91
92Revision 1.6 2000/06/28 12:19:18 morsch
93More consequent seperation of global input data services (AliMUONClusterInput singleton) and the
94cluster and hit reconstruction algorithms in AliMUONClusterFinderVS.
95AliMUONClusterFinderVS becomes the base class for clustering and hit reconstruction.
96It requires two cathode planes. Small modifications in the code will make it usable for
97one cathode plane and, hence, more general (for test beam data).
98AliMUONClusterFinder is now obsolete.
99
100Revision 1.5 2000/06/28 08:06:10 morsch
101Avoid global variables in AliMUONClusterFinderVS by seperating the input data for the fit from the
102algorithmic part of the class. Input data resides inside the AliMUONClusterInput singleton.
103It also naturally takes care of the TMinuit instance.
104
105Revision 1.4 2000/06/27 16:18:47 gosset
106Finally correct implementation of xm, ym, ixm, iym sizes
107when at least three local maxima on cathode 1 or on cathode 2
108
109Revision 1.3 2000/06/22 14:02:45 morsch
110Parameterised size of xm[], ym[], ixm[], iym[] correctly implemented (PH)
111Some HP scope problems corrected (PH)
112
113Revision 1.2 2000/06/15 07:58:48 morsch
114Code from MUON-dev joined
115
116Revision 1.1.2.3 2000/06/09 21:58:33 morsch
117Most coding rule violations corrected.
118
119Revision 1.1.2.2 2000/02/15 08:33:52 morsch
120Error in calculation of contribution map for double clusters (Split method) corrected (A.M.)
121Error in determination of track list for double cluster (FillCluster method) corrected (A.M.)
122Revised 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.
132
133*/
134
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"
142#include "AliRun.h"
143#include "AliMUON.h"
144
145#include <TTree.h>
146#include <TCanvas.h>
147#include <TH1.h>
148#include <TPad.h>
149#include <TGraph.h>
150#include <TPostScript.h>
151#include <TMinuit.h>
152#include <TF1.h>
153
154#include <stdio.h>
155#include <iostream.h>
156
157//_____________________________________________________________________
158// This function is minimized in the double-Mathieson fit
159void fcnS2(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag);
160void fcnS1(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag);
161void fcnCombiS1(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag);
162void fcnCombiS2(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag);
163
164ClassImp(AliMUONClusterFinderVS)
165
166 AliMUONClusterFinderVS::AliMUONClusterFinderVS()
167{
168// Default constructor
169 fInput=AliMUONClusterInput::Instance();
170 fHitMap[0] = 0;
171 fHitMap[1] = 0;
172 fTrack[0]=fTrack[1]=-1;
173 fDebugLevel = 0; // make silent default
174 fGhostChi2Cut = 1e6; // nothing done by default
175 fSeg[0] = 0;
176 fSeg[1] = 0;
177 for(Int_t i=0; i<100; i++) {
178 for (Int_t j=0; j<2; j++) {
179 fDig[i][j] = 0;
180 }
181 }
182}
183
184AliMUONClusterFinderVS::AliMUONClusterFinderVS(
185 const AliMUONClusterFinderVS & clusterFinder)
186{
187// Dummy copy Constructor
188 ;
189}
190
191void AliMUONClusterFinderVS::Decluster(AliMUONRawCluster *cluster)
192{
193// Decluster by local maxima
194 SplitByLocalMaxima(cluster);
195}
196
197void AliMUONClusterFinderVS::SplitByLocalMaxima(AliMUONRawCluster *c)
198{
199// Split complex cluster by local maxima
200 Int_t cath, i;
201
202 fInput->SetCluster(c);
203
204 fMul[0]=c->fMultiplicity[0];
205 fMul[1]=c->fMultiplicity[1];
206
207//
208// dump digit information into arrays
209//
210
211 Float_t qtot;
212
213 for (cath=0; cath<2; cath++) {
214 qtot=0;
215 for (i=0; i<fMul[cath]; i++)
216 {
217 // pointer to digit
218 fDig[i][cath]=fInput->Digit(cath, c->fIndexMap[i][cath]);
219 // pad coordinates
220 fIx[i][cath]= fDig[i][cath]->PadX();
221 fIy[i][cath]= fDig[i][cath]->PadY();
222 // pad charge
223 fQ[i][cath] = fDig[i][cath]->Signal();
224 // pad centre coordinates
225 fSeg[cath]->
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
229
230
231 FindLocalMaxima(c);
232
233//
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)
243
244 // One local maximum on cathodes 1 and 2 (X->cathode 2, Y->cathode 1)
245 if (fNLocal[0]==1 && fNLocal[1]==1) {
246 fXInit[0]=c->fX[1];
247 fYInit[0]=c->fY[0];
248 // One local maximum on cathode 1 (X,Y->cathode 1)
249 } else if (fNLocal[0]==1) {
250 fXInit[0]=c->fX[0];
251 fYInit[0]=c->fY[0];
252 // One local maximum on cathode 2 (X,Y->cathode 2)
253 } else {
254 fXInit[0]=c->fX[1];
255 fYInit[0]=c->fY[1];
256 }
257 if (fDebugLevel)
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);
264 oldchi2=chi2;
265 if (fDebugLevel)
266 fprintf(stderr," chi2 %f ",chi2);
267
268 c->fX[0]=fXFit[0];
269 c->fY[0]=fYFit[0];
270
271 c->fX[1]=fXFit[0];
272 c->fY[1]=fYFit[0];
273 c->fChi2[0]=chi2;
274 c->fChi2[1]=chi2;
275 // Force on anod
276 c->fX[0]=fSeg[0]->GetAnod(c->fX[0]);
277 c->fX[1]=fSeg[1]->GetAnod(c->fX[1]);
278
279// If reasonable chi^2 add result to the list of rawclusters
280 if (chi2 < 0.3) {
281 AddRawCluster(*c);
282// If not try combined double Mathieson Fit
283 } else {
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];
295 } else {
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];
300 }
301
302// Initial value for charge ratios
303 fQrInit[0]=0.5;
304 fQrInit[1]=0.5;
305 if (fDebugLevel)
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);
312
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
318 Split(c);
319 } else {
320 fprintf(stderr," Don't Split\n");
321 // Don't split
322 AddRawCluster(*c);
323 }
324 }
325
326// +++++++++++++++++++++++++++++++++++++++
327// (2) Two local maxima per cathode plane
328// +++++++++++++++++++++++++++++++++++++++
329 } else if (fNLocal[0]==2 && fNLocal[1]==2) {
330//
331// Let's look for ghosts first
332
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;
337//
338// Form the 2x2 combinations
339// 0-0, 0-1, 1-0, 1-1
340 ico=0;
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];
347
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];
352 ico++;
353 }
354 }
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
359
360// Analyse the combinations and keep those that are possible !
361// For each combination check consistency in x and y
362 Int_t iacc;
363 Bool_t accepted[4];
364 Float_t dr[4] = {1.e4, 1.e4, 1.e4, 1.e4};
365 iacc=0;
366
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;
369 Float_t eps = 1.e-5;
370//
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]);
381 if (fDebugLevel>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)) {
384 // consistent
385 accepted[ico]=kTRUE;
386 dr[ico] = TMath::Sqrt(dx*dx+dy*dy);
387 iacc++;
388 } else {
389 // reject
390 accepted[ico]=kFALSE;
391 }
392 }
393 printf("\n iacc= %d:\n", iacc);
394 if (iacc == 3) {
395 if (accepted[0] && accepted[1]) {
396 if (dr[0] >= dr[1]) {
397 accepted[0]=kFALSE;
398 } else {
399 accepted[1]=kFALSE;
400 }
401 }
402
403 if (accepted[2] && accepted[3]) {
404 if (dr[2] >= dr[3]) {
405 accepted[2]=kFALSE;
406 } else {
407 accepted[3]=kFALSE;
408 }
409 }
410/*
411// eliminate one candidate
412 Float_t drmax = 0;
413 Int_t icobad = -1;
414
415 for (ico=0; ico<4; ico++) {
416 if (accepted[ico] && dr[ico] > drmax) {
417 icobad = ico;
418 drmax = dr[ico];
419 }
420 }
421
422 accepted[icobad] = kFALSE;
423*/
424 iacc = 2;
425 }
426
427
428 printf("\n iacc= %d:\n", iacc);
429 if (fDebugLevel) {
430 if (iacc==2) {
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");
436 }
437 }
438
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]);
444
445// ******* iacc = 0 *******
446// No combinations found between the 2 cathodes
447// We keep the center of gravity of the cluster
448 if (iacc==0) {
449 AddRawCluster(*c);
450 }
451
452// ******* iacc = 1 *******
453// Only one combination found between the 2 cathodes
454 if (iacc==1) {
455// Initial values for the 2 maxima (x,y)
456
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
459 if (accepted[0]){
460 fprintf(stderr,"ico=0\n");
461 fXInit[0]=xm[0][1];
462 fYInit[0]=ym[0][0];
463 fXInit[1]=xm[3][0];
464 fYInit[1]=ym[3][0];
465 } else if (accepted[1]){
466 fprintf(stderr,"ico=1\n");
467 fXInit[0]=xm[1][1];
468 fYInit[0]=ym[1][0];
469 fXInit[1]=xm[2][0];
470 fYInit[1]=ym[2][0];
471 } else if (accepted[2]){
472 fprintf(stderr,"ico=2\n");
473 fXInit[0]=xm[2][1];
474 fYInit[0]=ym[2][0];
475 fXInit[1]=xm[1][0];
476 fYInit[1]=ym[1][0];
477 } else if (accepted[3]){
478 fprintf(stderr,"ico=3\n");
479 fXInit[0]=xm[3][1];
480 fYInit[0]=ym[3][0];
481 fXInit[1]=xm[0][0];
482 fYInit[1]=ym[0][0];
483 }
484 if (fDebugLevel)
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);
491 if (fDebugLevel)
492 fprintf(stderr," chi2 %f\n",chi2);
493
494// If reasonable chi^2 add result to the list of rawclusters
495 if (chi2<10) {
496 Split(c);
497
498 } else {
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
501 if (accepted[0]){
502 fprintf(stderr,"ico=0\n");
503 fXInit[0]=xm[0][1];
504 fYInit[0]=ym[0][0];
505 fXInit[1]=xm[3][1];
506 fYInit[1]=ym[3][1];
507 } else if (accepted[1]){
508 fprintf(stderr,"ico=1\n");
509 fXInit[0]=xm[1][1];
510 fYInit[0]=ym[1][0];
511 fXInit[1]=xm[2][1];
512 fYInit[1]=ym[2][1];
513 } else if (accepted[2]){
514 fprintf(stderr,"ico=2\n");
515 fXInit[0]=xm[2][1];
516 fYInit[0]=ym[2][0];
517 fXInit[1]=xm[1][1];
518 fYInit[1]=ym[1][1];
519 } else if (accepted[3]){
520 fprintf(stderr,"ico=3\n");
521 fXInit[0]=xm[3][1];
522 fYInit[0]=ym[3][0];
523 fXInit[1]=xm[0][1];
524 fYInit[1]=ym[0][1];
525 }
526 if (fDebugLevel)
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);
533 if (fDebugLevel)
534 fprintf(stderr," chi2 %f\n",chi2);
535
536// If reasonable chi^2 add result to the list of rawclusters
537 if (chi2<10) {
538 Split(c);
539 } else {
540//We keep only the combination found (X->cathode 2, Y->cathode 1)
541 for (Int_t ico=0; ico<2; ico++) {
542 if (accepted[ico]) {
543 AliMUONRawCluster cnew;
544 Int_t cath;
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;
549
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]);
554 }
555 fprintf(stderr,"\nRawCluster %d cath %d\n",ico,cath);
556 fprintf(stderr,"mult_av %d\n",c->fMultiplicity[cath]);
557 FillCluster(&cnew,cath);
558 }
559 cnew.fClusterType=cnew.PhysicsContribution();
560 AddRawCluster(cnew);
561 fNPeaks++;
562 }
563 }
564 }
565 }
566 }
567
568// ******* iacc = 2 *******
569// Two combinations found between the 2 cathodes
570 if (iacc==2) {
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");
574
575// Have a try !! with that
576 if (accepted[0]&&accepted[3]) {
577 fXInit[0]=xm[0][1];
578 fYInit[0]=ym[0][0];
579 fXInit[1]=xm[1][1];
580 fYInit[1]=ym[1][0];
581 } else {
582 fXInit[0]=xm[2][1];
583 fYInit[0]=ym[2][0];
584 fXInit[1]=xm[3][1];
585 fYInit[1]=ym[3][0];
586 }
587 if (fDebugLevel)
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);
594 Split(c);
595
596 } else {
597// No ghosts ! No Problems ! - Perform one fit only !
598 if (accepted[0]&&accepted[3]) {
599 fXInit[0]=xm[0][1];
600 fYInit[0]=ym[0][0];
601 fXInit[1]=xm[3][1];
602 fYInit[1]=ym[3][0];
603 } else {
604 fXInit[0]=xm[1][1];
605 fYInit[0]=ym[1][0];
606 fXInit[1]=xm[2][1];
607 fYInit[1]=ym[2][0];
608 }
609 if (fDebugLevel)
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);
616 if (fDebugLevel)
617 fprintf(stderr," chi2 %f\n",chi2);
618 Split(c);
619 }
620
621// ******* iacc = 4 *******
622// Four combinations found between the 2 cathodes
623// Ghost !!
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
628 fXInit[0]=xm[0][1];
629 fYInit[0]=ym[0][0];
630 fXInit[1]=xm[3][1];
631 fYInit[1]=ym[3][0];
632 if (fDebugLevel)
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);
639 if (fDebugLevel)
640 fprintf(stderr," chi2 %f\n",chi2);
641 // store results of fit and postpone decision
642 Double_t sXFit[2],sYFit[2],sQrFit[2];
643 Float_t sChi2[2];
644 for (Int_t i=0;i<2;i++) {
645 sXFit[i]=fXFit[i];
646 sYFit[i]=fYFit[i];
647 sQrFit[i]=fQrFit[i];
648 sChi2[i]=fChi2[i];
649 }
650 fXInit[0]=xm[1][1];
651 fYInit[0]=ym[1][0];
652 fXInit[1]=xm[2][1];
653 fYInit[1]=ym[2][0];
654 if (fDebugLevel)
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);
661 if (fDebugLevel)
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;
666
667 chi2f = (TMath::Log(fInput->TotalCharge(0)*fQrFit[0]
668 / (fInput->TotalCharge(1)*fQrFit[1]) )
669 / fInput->Response()->ChargeCorrel() );
670 chi2f *=chi2f;
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;
675
676 chi2s = (TMath::Log(fInput->TotalCharge(0)*sQrFit[0]
677 / (fInput->TotalCharge(1)*sQrFit[1]) )
678 / fInput->Response()->ChargeCorrel() );
679 chi2s *=chi2s;
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;
684
685 // usefull to store the charge matching chi2 in the cluster
686 // fChi2[0]=sChi2[1]=chi2f;
687 // fChi2[1]=sChi2[0]=chi2s;
688
689 if (chi2f<=fGhostChi2Cut && chi2s<=fGhostChi2Cut)
690 c->fGhost=1;
691 if (chi2f>fGhostChi2Cut && chi2s>fGhostChi2Cut) {
692 // we keep the ghost
693 c->fGhost=2;
694 chi2s=-1;
695 chi2f=-1;
696 }
697 if (chi2f<=fGhostChi2Cut)
698 Split(c);
699 if (chi2s<=fGhostChi2Cut) {
700 // retreive saved values
701 for (Int_t i=0;i<2;i++) {
702 fXFit[i]=sXFit[i];
703 fYFit[i]=sYFit[i];
704 fQrFit[i]=sQrFit[i];
705 fChi2[i]=sChi2[i];
706 }
707 Split(c);
708 }
709 c->fGhost=0;
710 }
711
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// ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
716//
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;
721//
722// Form the 2x2 combinations
723// 0-0, 0-1, 1-0, 1-1
724 ico=0;
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];
730
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];
735 ico++;
736 }
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
739
740// Analyse the combinations and keep those that are possible !
741// For each combination check consistency in x and y
742 Int_t iacc;
743 Bool_t accepted[4];
744 iacc=0;
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;
747
748 Float_t eps = 1.e-5;
749
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]);
758 if (fDebugLevel>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)) {
761 // consistent
762 accepted[ico]=kTRUE;
763 iacc++;
764 } else {
765 // reject
766 accepted[ico]=kFALSE;
767 }
768 }
769
770 Float_t chi21 = 100;
771 Float_t chi22 = 100;
772 Float_t chi23 = 100;
773
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];
778
779 if (accepted[0] && accepted[1]) {
780
781 fXInit[0]=0.5*(xm[0][1]+xm[0][0]);
782 fYInit[0]=ym[0][0];
783 fXInit[1]=0.5*(xm[0][1]+xm[1][0]);
784 fYInit[1]=ym[1][0];
785 fQrInit[0]=0.5;
786 fQrInit[1]=0.5;
787 chi23=CombiDoubleMathiesonFit(c);
788 if (chi23<10) {
789 Split(c);
790 Float_t yst;
791 yst = fYFit[0];
792 fYFit[0] = fYFit[1];
793 fYFit[1] = yst;
794 Split(c);
795 }
796 } else if (accepted[0]) {
797 fXInit[0]=xm[0][1];
798 fYInit[0]=ym[0][0];
799 fXInit[1]=xm[1][0];
800 fYInit[1]=ym[1][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);
806 if (fDebugLevel)
807 fprintf(stderr," chi2 %f\n",chi21);
808 if (chi21<10) Split(c);
809 } else if (accepted[1]) {
810 fXInit[0]=xm[1][1];
811 fYInit[0]=ym[1][0];
812 fXInit[1]=xm[0][0];
813 fYInit[1]=ym[0][0];
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);
819 if (fDebugLevel)
820 fprintf(stderr," chi2 %f\n",chi22);
821 if (chi22<10) Split(c);
822 }
823
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++) {
827 if (accepted[ico]) {
828 AliMUONRawCluster cnew;
829 Int_t cath;
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]);
838 }
839 fprintf(stderr,"\nRawCluster %d cath %d\n",ico,cath);
840 fprintf(stderr,"mult_av %d\n",c->fMultiplicity[cath]);
841 FillCluster(&cnew,cath);
842 }
843 cnew.fClusterType=cnew.PhysicsContribution();
844 AddRawCluster(cnew);
845 fNPeaks++;
846 }
847 }
848 }
849
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;
858//
859// Form the 2x2 combinations
860// 0-0, 0-1, 1-0, 1-1
861 ico=0;
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];
867
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];
872 ico++;
873 }
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
876
877// Analyse the combinations and keep those that are possible !
878// For each combination check consistency in x and y
879 Int_t iacc;
880 Bool_t accepted[4];
881 iacc=0;
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;
884 Float_t eps = 1.e-5;
885
886
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]);
895 if (fDebugLevel>0)
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)) {
898 // consistent
899 accepted[ico]=kTRUE;
900 fprintf(stderr,"ico %d\n",ico);
901 iacc++;
902 } else {
903 // reject
904 accepted[ico]=kFALSE;
905 }
906 }
907
908 Float_t chi21 = 100;
909 Float_t chi22 = 100;
910 Float_t chi23 = 100;
911
912 fQrInit[1]=Float_t(fQ[fIndLocal[0][1]][1])/
913 Float_t(fQ[fIndLocal[0][1]][1]+fQ[fIndLocal[1][1]][1]);
914
915 fQrInit[0]=fQrInit[1];
916
917
918 if (accepted[0] && accepted[1]) {
919 fXInit[0]=xm[0][1];
920 fYInit[0]=0.5*(ym[0][0]+ym[0][1]);
921 fXInit[1]=xm[1][1];
922 fYInit[1]=0.5*(ym[0][0]+ym[1][1]);
923 fQrInit[0]=0.5;
924 fQrInit[1]=0.5;
925 chi23=CombiDoubleMathiesonFit(c);
926 if (chi23<10) {
927 Split(c);
928 Float_t yst;
929 yst = fYFit[0];
930 fYFit[0] = fYFit[1];
931 fYFit[1] = yst;
932 Split(c);
933 }
934 } else if (accepted[0]) {
935 fXInit[0]=xm[0][0];
936 fYInit[0]=ym[0][1];
937 fXInit[1]=xm[1][1];
938 fYInit[1]=ym[1][1];
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);
944 if (fDebugLevel)
945 fprintf(stderr," chi2 %f\n",chi21);
946 if (chi21<10) Split(c);
947 } else if (accepted[1]) {
948 fXInit[0]=xm[1][0];
949 fYInit[0]=ym[1][1];
950 fXInit[1]=xm[0][1];
951 fYInit[1]=ym[0][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);
957 if (fDebugLevel)
958 fprintf(stderr," chi2 %f\n",chi22);
959 if (chi22<10) Split(c);
960 }
961
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++) {
965 if (accepted[ico]) {
966 AliMUONRawCluster cnew;
967 Int_t cath;
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]);
976 }
977 fprintf(stderr,"\nRawCluster %d cath %d\n",ico,cath);
978 fprintf(stderr,"mult_av %d\n",c->fMultiplicity[cath]);
979 FillCluster(&cnew,cath);
980 }
981 cnew.fClusterType=cnew.PhysicsContribution();
982 AddRawCluster(cnew);
983 fNPeaks++;
984 }
985 }
986 }
987
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];
993 Int_t ii;
994
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];
1003
1004 Int_t isec, ico;
1005 Float_t dpx, dpy, dx, dy;
1006
1007 ico=0;
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];
1014
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];
1019 ico++;
1020 }
1021 }
1022
1023 Int_t nIco = ico;
1024 if (fDebugLevel)
1025 fprintf(stderr,"nIco %d\n",nIco);
1026 for (ico=0; ico<nIco; ico++) {
1027 if (fDebugLevel)
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]);
1035 if (fDebugLevel) {
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]);
1038 }
1039 if ((dx <= dpx) && (dy <= dpy)) {
1040 if (fDebugLevel)
1041 fprintf(stderr,"ok\n");
1042 Int_t cath;
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]);
1052 }
1053 FillCluster(&cnew,cath);
1054 }
1055 cnew.fClusterType=cnew.PhysicsContribution();
1056 AddRawCluster(cnew);
1057 fNPeaks++;
1058 }
1059 }
1060 delete [] xm;
1061 delete [] ym;
1062 delete [] ixm;
1063 delete [] iym;
1064 }
1065}
1066
1067void AliMUONClusterFinderVS::FindLocalMaxima(AliMUONRawCluster* c)
1068{
1069// Find all local maxima of a cluster
1070 if (fDebugLevel)
1071 printf("\n Find Local maxima !");
1072
1073 AliMUONDigit* digt;
1074
1075 Int_t cath, cath1; // loops over cathodes
1076 Int_t i; // loops over digits
1077 Int_t j; // loops over cathodes
1078//
1079// Find local maxima
1080//
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;
1087 }
1088// number of next neighbours and arrays to store them
1089 Int_t nn;
1090 Int_t x[10], y[10];
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;
1109 break;
1110//
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])
1117 {
1118 isLocal[i][cath]=kFALSE;
1119 }
1120 } // loop over local maxima
1121 } // are there already local maxima
1122 } // same charge ?
1123 } // loop over next neighbours
1124 if (isLocal[i][cath]) {
1125 fIndLocal[fNLocal[cath]][cath]=i;
1126 fNLocal[cath]++;
1127 }
1128 } // loop over all digits
1129 } // loop over cathodes
1130
1131 if (fDebugLevel) {
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]);
1136 }
1137 Int_t ix, iy, isec;
1138 Float_t dpx, dpy;
1139
1140
1141 if (fNLocal[1]==2 && (fNLocal[0]==1 || fNLocal[0]==0)) {
1142 Int_t iback=fNLocal[0];
1143
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
1146//
1147// Loop over cluster digits
1148 cath=0;
1149 cath1=1;
1150
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.))
1159 continue;
1160
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
1165 Int_t iNN=0;
1166
1167 for (fSeg[cath]
1168 ->FirstPad(fX[i][cath], fY[i][cath], fZPlane, 0., dpy);
1169 fSeg[cath]
1170 ->MorePads();
1171 fSeg[cath]
1172 ->NextPad())
1173 {
1174 ix = fSeg[cath]->Ix();
1175 iy = fSeg[cath]->Iy();
1176 // skip the current pad
1177 if (iy == fIy[i][cath]) continue;
1178
1179 if (fHitMap[cath]->TestHit(ix, iy)!=kEmpty) {
1180 iNN++;
1181 digt=(AliMUONDigit*) fHitMap[cath]->GetHit(ix,iy);
1182 if (digt->Signal() > fQ[i][cath]) isLocal[i][cath]=kFALSE;
1183 }
1184 } // Loop over pad neighbours in y
1185 if (isLocal[i][cath] && iNN>0) {
1186 fIndLocal[fNLocal[cath]][cath]=i;
1187 fNLocal[cath]++;
1188 }
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
1192 if (fDebugLevel) {
1193 fprintf(stderr,
1194 "\n New search gives %d local maxima for cathode 1 \n",
1195 fNLocal[0]);
1196 fprintf(stderr,
1197 " %d local maxima for cathode 2 \n",
1198 fNLocal[1]);
1199 }
1200 if (fNLocal[cath]>2) {
1201 fNLocal[cath]=iback;
1202 }
1203
1204 } // 1,2 local maxima
1205
1206 if (fNLocal[0]==2 && (fNLocal[1]==1 || fNLocal[1]==0)) {
1207 Int_t iback=fNLocal[1];
1208
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
1211 cath=1;
1212 Int_t cath1 = 0;
1213 Float_t eps = 1.e-5;
1214
1215//
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))
1225 continue;
1226
1227//
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
1231
1232// iNN counts the number of neighbours with signal, it should be 1 or 2
1233 Int_t iNN=0;
1234 for (fSeg[cath]
1235 ->FirstPad(fX[i][cath], fY[i][cath], fZPlane, dpx, 0.);
1236 fSeg[cath]
1237 ->MorePads();
1238 fSeg[cath]
1239 ->NextPad())
1240 {
1241
1242 ix = fSeg[cath]->Ix();
1243 iy = fSeg[cath]->Iy();
1244
1245 // skip the current pad
1246 if (ix == fIx[i][cath]) continue;
1247
1248 if (fHitMap[cath]->TestHit(ix, iy)!=kEmpty) {
1249 iNN++;
1250 digt=(AliMUONDigit*) fHitMap[cath]->GetHit(ix,iy);
1251 if (digt->Signal() > fQ[i][cath]) isLocal[i][cath]=kFALSE;
1252 }
1253 } // Loop over pad neighbours in x
1254 if (isLocal[i][cath] && iNN>0) {
1255 fIndLocal[fNLocal[cath]][cath]=i;
1256 fNLocal[cath]++;
1257 }
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
1261 if (fDebugLevel) {
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]);
1265 }
1266 if (fNLocal[cath]>2) {
1267 fNLocal[cath]=iback;
1268 }
1269 } // 2,1 local maxima
1270}
1271
1272
1273void AliMUONClusterFinderVS::FillCluster(AliMUONRawCluster* c, Int_t flag, Int_t cath)
1274{
1275//
1276// Completes cluster information starting from list of digits
1277//
1278 AliMUONDigit* dig;
1279 Float_t x, y, z;
1280 Int_t ix, iy;
1281
1282 if (cath==1) {
1283 c->fPeakSignal[cath]=c->fPeakSignal[0];
1284 } else {
1285 c->fPeakSignal[cath]=0;
1286 }
1287
1288
1289 if (flag) {
1290 c->fX[cath]=0;
1291 c->fY[cath]=0;
1292 c->fQ[cath]=0;
1293 }
1294
1295 if (fDebugLevel)
1296 fprintf(stderr,"\n fPeakSignal %d\n",c->fPeakSignal[cath]);
1297 for (Int_t i=0; i<c->fMultiplicity[cath]; i++)
1298 {
1299 dig= fInput->Digit(cath,c->fIndexMap[i][cath]);
1300 ix=dig->PadX()+c->fOffsetMap[i][cath];
1301 iy=dig->PadY();
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;
1310//
1311//
1312 if (fDebugLevel>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]);
1321 }
1322//
1323 if (flag) {
1324 fSeg[cath]->GetPadC(ix, iy, x, y, z);
1325 c->fX[cath] += q*x;
1326 c->fY[cath] += q*y;
1327 c->fQ[cath] += q;
1328 }
1329 } // loop over digits
1330 if (fDebugLevel)
1331 fprintf(stderr," fin du cluster c\n");
1332
1333
1334 if (flag) {
1335 c->fX[cath]/=c->fQ[cath];
1336// Force on anod
1337 c->fX[cath]=fSeg[cath]->GetAnod(c->fX[cath]);
1338 c->fY[cath]/=c->fQ[cath];
1339//
1340// apply correction to the coordinate along the anode wire
1341//
1342 x=c->fX[cath];
1343 y=c->fY[cath];
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);
1348
1349 if (cogCorr) {
1350 Float_t yOnPad=(c->fY[cath]-y)/fSeg[cath]->Dpy(isec);
1351 c->fY[cath]=c->fY[cath]-cogCorr->Eval(yOnPad, 0, 0);
1352 }
1353 }
1354}
1355
1356void AliMUONClusterFinderVS::FillCluster(AliMUONRawCluster* c, Int_t cath)
1357{
1358//
1359// Completes cluster information starting from list of digits
1360//
1361 static Float_t dr0;
1362
1363 AliMUONDigit* dig;
1364
1365 if (cath==0) {
1366 dr0 = 10000;
1367 }
1368
1369 Float_t xpad, ypad, zpad;
1370 Float_t dx, dy, dr;
1371
1372 for (Int_t i=0; i<c->fMultiplicity[cath]; i++)
1373 {
1374 dig = fInput->Digit(cath,c->fIndexMap[i][cath]);
1375 fSeg[cath]->
1376 GetPadC(dig->PadX(),dig->PadY(),xpad,ypad, zpad);
1377 if (fDebugLevel)
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);
1382
1383 if (dr < dr0) {
1384 dr0 = dr;
1385 if (fDebugLevel)
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);
1397 if (fDebugLevel)
1398 fprintf(stderr," c->fTracks[0] %d c->fTracks[1] %d\n",dig->Hit(),
1399 dig->Track(0));
1400 }
1401//
1402 } // loop over digits
1403
1404// apply correction to the coordinate along the anode wire
1405// Force on anod
1406 c->fX[cath]=fSeg[cath]->GetAnod(c->fX[cath]);
1407}
1408
1409void AliMUONClusterFinderVS::FindCluster(Int_t i, Int_t j, Int_t cath, AliMUONRawCluster &c){
1410
1411
1412//
1413// Find a super cluster on both cathodes
1414//
1415//
1416// Add i,j as element of the cluster
1417//
1418
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();
1424
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);
1430 }
1431
1432//
1433// Make sure that list of digits is ordered
1434//
1435 Int_t mu=c.fMultiplicity[cath];
1436 c.fIndexMap[mu][cath]=idx;
1437
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;
1443
1444
1445 if (mu > 0) {
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();
1451
1452 if (q>ql || (q==ql && theX > ix && theY < iy)) {
1453 c.fIndexMap[ind][cath]=idx;
1454 c.fIndexMap[ind+1][cath]=ist;
1455 } else {
1456
1457 break;
1458 }
1459 }
1460 }
1461
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;
1466 }
1467
1468// Prepare center of gravity calculation
1469 Float_t x, y, z;
1470 fSeg[cath]->GetPadC(i, j, x, y, z);
1471
1472 c.fX[cath] += q*x;
1473 c.fY[cath] += q*y;
1474 c.fQ[cath] += q;
1475//
1476// Flag hit as "taken"
1477 fHitMap[cath]->FlagHit(i,j);
1478//
1479// Now look recursively for all neighbours and pad hit on opposite cathode
1480//
1481// Loop over neighbours
1482 Int_t ix,iy;
1483 ix=iy=0;
1484 Int_t nn;
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++) {
1488 ix=xList[in];
1489 iy=yList[in];
1490
1491 if (fHitMap[cath]->TestHit(ix,iy)==kUnused) {
1492 if (fDebugLevel>1)
1493 printf("\n Neighbours %d %d %d", cath, ix, iy);
1494 FindCluster(ix, iy, cath, c);
1495 }
1496
1497 }
1498 Int_t nOpp=0;
1499 Int_t iXopp[50], iYopp[50];
1500
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);
1504 Int_t iop;
1505 Float_t dx, dy;
1506
1507 if (cath==0) {
1508 iop = 1;
1509 dx = (fSeg[cath]->Dpx(isec))/2.;
1510 dy = 0.;
1511 } else {
1512 iop = 0;
1513 dx = 0.;
1514 dy = (fSeg[cath]->Dpy(isec))/2;
1515 }
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())
1520 {
1521
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){
1526 iXopp[nOpp]=ix;
1527 iYopp[nOpp++]=iy;
1528 if (fDebugLevel > 1)
1529 printf("\n Opposite %d %d %d", iop, ix, iy);
1530 }
1531
1532 } // Loop over pad neighbours
1533// This had to go outside the loop since recursive calls inside the iterator are not possible
1534//
1535 Int_t jopp;
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);
1539 }
1540}
1541
1542//_____________________________________________________________________________
1543
1544void AliMUONClusterFinderVS::FindRawClusters()
1545{
1546 //
1547 // MUON cluster finder from digits -- finds neighbours on both cathodes and
1548 // fills the tree with raw clusters
1549 //
1550
1551// Return if no input datad available
1552 if (!fInput->NDigits(0) && !fInput->NDigits(1)) return;
1553
1554 fSeg[0] = fInput->Segmentation(0);
1555 fSeg[1] = fInput->Segmentation(1);
1556
1557 fHitMap[0] = new AliMUONHitMapA1(fSeg[0], fInput->Digits(0));
1558 fHitMap[1] = new AliMUONHitMapA1(fSeg[1], fInput->Digits(1));
1559
1560
1561 AliMUONDigit *dig;
1562
1563 Int_t ndig, cath;
1564 Int_t nskip=0;
1565 Int_t ncls=0;
1566 fHitMap[0]->FillHits();
1567 fHitMap[1]->FillHits();
1568//
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) {
1576 nskip++;
1577 continue;
1578 }
1579 if (fDebugLevel)
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
1589 c.fNcluster[0]=-1;
1590 Float_t xcu, ycu;
1591 fSeg[cath]->GetPadC(i,j,xcu, ycu, fZPlane);
1592 fSector= fSeg[cath]->Sector(i,j)/100;
1593 if (fDebugLevel)
1594 printf("\n New Seed %d %d ", i,j);
1595
1596 FindCluster(i,j,cath,c);
1597// ^^^^^^^^^^^^^^^^^^^^^^^^
1598 // center of gravity
1599 c.fX[0] /= c.fQ[0];
1600// Force on anod
1601 c.fX[0]=fSeg[0]->GetAnod(c.fX[0]);
1602 c.fY[0] /= c.fQ[0];
1603 c.fX[1] /= c.fQ[1];
1604// Force on anod
1605 c.fX[1]=fSeg[0]->GetAnod(c.fX[1]);
1606 c.fY[1] /= c.fQ[1];
1607
1608 c.fZ[0] = fZPlane;
1609 c.fZ[1] = fZPlane;
1610
1611 if (fDebugLevel) {
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]);
1616 }
1617// Analyse cluster and decluster if necessary
1618//
1619 ncls++;
1620 c.fNcluster[1]=fNRawClusters;
1621 c.fClusterType=c.PhysicsContribution();
1622
1623 fNPeaks=0;
1624//
1625//
1626 Decluster(&c);
1627//
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
1632
1633 { // begin local scope
1634 for (int k=0;k<c.fMultiplicity[1];k++) c.fIndexMap[k][1]=0;
1635 } // end local scope
1636
1637 c.fMultiplicity[0]=c.fMultiplicity[0]=0;
1638
1639
1640 } // end loop ndig
1641 } // end loop cathodes
1642 delete fHitMap[0];
1643 delete fHitMap[1];
1644}
1645
1646Float_t AliMUONClusterFinderVS::SingleMathiesonFit(AliMUONRawCluster *c, Int_t cath)
1647{
1648// Performs a single Mathieson fit on one cathode
1649//
1650 AliMUONClusterInput& clusterInput = *(AliMUONClusterInput::Instance());
1651
1652 clusterInput.Fitter()->SetFCN(fcnS1);
1653 clusterInput.Fitter()->mninit(2,10,7);
1654 Double_t arglist[20];
1655 Int_t ierflag=0;
1656 arglist[0]=1;
1657// Set starting values
1658 static Double_t vstart[2];
1659 vstart[0]=c->fX[1];
1660 vstart[1]=c->fY[0];
1661
1662
1663// lower and upper limits
1664 static Double_t lower[2], upper[2];
1665 Int_t ix,iy;
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;
1670
1671 upper[0]=lower[0]+fSeg[cath]->Dpx(isec);
1672 upper[1]=lower[1]+fSeg[cath]->Dpy(isec);
1673
1674// step sizes
1675 static Double_t step[2]={0.0005, 0.0005};
1676
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);
1681 if (fDebugLevel==0)
1682 clusterInput.Fitter()->mnexcm("SET NOW", arglist, 0, ierflag);
1683 clusterInput.Fitter()->mnexcm("SET OUT", arglist, 0, ierflag);
1684 arglist[0]= -1;
1685 arglist[1]= 0;
1686
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;
1692
1693 clusterInput.Fitter()->mnstat(fmin, fedm, errdef, npari, nparx, istat);
1694 fFitStat=istat;
1695
1696// Print results
1697// Get fitted parameters
1698 Double_t xrec, yrec;
1699 TString chname;
1700 Double_t epxz, b1, b2;
1701 Int_t ierflg;
1702 clusterInput.Fitter()->mnpout(0, chname, xrec, epxz, b1, b2, ierflg);
1703 clusterInput.Fitter()->mnpout(1, chname, yrec, epxz, b1, b2, ierflg);
1704 fXFit[cath]=xrec;
1705 fYFit[cath]=yrec;
1706 return fmin;
1707}
1708
1709Float_t AliMUONClusterFinderVS::CombiSingleMathiesonFit(AliMUONRawCluster *c)
1710{
1711// Perform combined Mathieson fit on both cathode planes
1712//
1713 AliMUONClusterInput& clusterInput = *(AliMUONClusterInput::Instance());
1714 clusterInput.Fitter()->SetFCN(fcnCombiS1);
1715 clusterInput.Fitter()->mninit(2,10,7);
1716 Double_t arglist[20];
1717 Int_t ierflag=0;
1718 arglist[0]=1;
1719 static Double_t vstart[2];
1720 vstart[0]=fXInit[0];
1721 vstart[1]=fYInit[0];
1722
1723
1724// lower and upper limits
1725 static Float_t lower[2], upper[2];
1726 Int_t ix,iy,isec;
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);
1733
1734 Int_t icount;
1735 Float_t xdum, ydum, zdum;
1736
1737// Find save upper and lower limits
1738
1739 icount = 0;
1740
1741 for (fSeg[1]->FirstPad(fXInit[0], fYInit[0], fZPlane, dpx, 0.);
1742 fSeg[1]->MorePads(); fSeg[1]->NextPad())
1743 {
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];
1747 icount++;
1748 }
1749
1750 if (lower[0]>upper[0]) {xdum=lower[0]; lower[0]=upper[0]; upper[0]=xdum;}
1751
1752 icount=0;
1753 if (fDebugLevel)
1754 printf("\n single y %f %f", fXInit[0], fYInit[0]);
1755
1756 for (fSeg[0]->FirstPad(fXInit[0], fYInit[0], fZPlane, 0., dpy);
1757 fSeg[0]->MorePads(); fSeg[0]->NextPad())
1758 {
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];
1762 icount++;
1763 if (fDebugLevel)
1764 printf("\n upper lower %d %f %f", icount, upper[1], lower[1]);
1765 }
1766
1767 if (lower[1]>upper[1]) {xdum=lower[1]; lower[1]=upper[1]; upper[1]=xdum;}
1768
1769// step sizes
1770 static Double_t step[2]={0.00001, 0.0001};
1771
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);
1776 if (fDebugLevel==0)
1777 clusterInput.Fitter()->mnexcm("SET NOW", arglist, 0, ierflag);
1778 clusterInput.Fitter()->mnexcm("SET OUT", arglist, 0, ierflag);
1779 arglist[0]= -1;
1780 arglist[1]= 0;
1781
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;
1787
1788 clusterInput.Fitter()->mnstat(fmin, fedm, errdef, npari, nparx, istat);
1789 fFitStat=istat;
1790
1791// Print results
1792// Get fitted parameters
1793 Double_t xrec, yrec;
1794 TString chname;
1795 Double_t epxz, b1, b2;
1796 Int_t ierflg;
1797 clusterInput.Fitter()->mnpout(0, chname, xrec, epxz, b1, b2, ierflg);
1798 clusterInput.Fitter()->mnpout(1, chname, yrec, epxz, b1, b2, ierflg);
1799 fXFit[0]=xrec;
1800 fYFit[0]=yrec;
1801 return fmin;
1802}
1803
1804Bool_t AliMUONClusterFinderVS::DoubleMathiesonFit(AliMUONRawCluster *c, Int_t cath)
1805{
1806// Performs a double Mathieson fit on one cathode
1807//
1808
1809//
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];
1815 Int_t ierflag=0;
1816 arglist[0]=1;
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);
1830
1831 upper[0]=lower[0]+2.*fSeg[cath]->Dpx(isec);
1832 upper[1]=lower[1]+2.*fSeg[cath]->Dpy(isec);
1833
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;
1837
1838 upper[2]=lower[2]+fSeg[cath]->Dpx(isec);
1839 upper[3]=lower[3]+fSeg[cath]->Dpy(isec);
1840
1841 lower[4]=0.;
1842 upper[4]=1.;
1843// step sizes
1844 static Double_t step[5]={0.0005, 0.0005, 0.0005, 0.0005, 0.0001};
1845
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);
1853 if (fDebugLevel==0)
1854 clusterInput.Fitter()->mnexcm("SET NOW", arglist, 0, ierflag);
1855 clusterInput.Fitter()->mnexcm("SET OUT", arglist, 0, ierflag);
1856 arglist[0]= -1;
1857 arglist[1]= 0;
1858
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;
1864 TString chname;
1865 Double_t epxz, b1, b2;
1866 Int_t ierflg;
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);
1872
1873 Double_t fmin, fedm, errdef;
1874 Int_t npari, nparx, istat;
1875
1876 clusterInput.Fitter()->mnstat(fmin, fedm, errdef, npari, nparx, istat);
1877 fFitStat=istat;
1878 return kTRUE;
1879}
1880
1881Float_t AliMUONClusterFinderVS::CombiDoubleMathiesonFit(AliMUONRawCluster *c)
1882{
1883//
1884// Perform combined double Mathieson fit on both cathode planes
1885//
1886 AliMUONClusterInput& clusterInput = *(AliMUONClusterInput::Instance());
1887 clusterInput.Fitter()->SetFCN(fcnCombiS2);
1888 clusterInput.Fitter()->mninit(6,10,7);
1889 Double_t arglist[20];
1890 Int_t ierflag=0;
1891 arglist[0]=1;
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];
1902 Int_t ix,iy,isec;
1903 Float_t dpx, dpy;
1904
1905 fSeg[1]->GetPadI(fXInit[0], fYInit[0], fZPlane, ix, iy);
1906 isec=fSeg[1]->Sector(ix, iy);
1907 dpx=fSeg[1]->Dpx(isec);
1908
1909 fSeg[0]->GetPadI(fXInit[0], fYInit[0], fZPlane, ix, iy);
1910 isec=fSeg[0]->Sector(ix, iy);
1911 dpy=fSeg[0]->Dpy(isec);
1912
1913
1914 Int_t icount;
1915 Float_t xdum, ydum, zdum;
1916 if (fDebugLevel)
1917 printf("\n Cluster Finder: %f %f %f %f ", fXInit[0], fXInit[1],fYInit[0], fYInit[1] );
1918
1919// Find save upper and lower limits
1920 icount = 0;
1921
1922 for (fSeg[1]->FirstPad(fXInit[0], fYInit[0], fZPlane, dpx, 0.);
1923 fSeg[1]->MorePads(); fSeg[1]->NextPad())
1924 {
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];
1929 icount++;
1930 }
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]);
1933
1934
1935 icount=0;
1936
1937 for (fSeg[0]->FirstPad(fXInit[0], fYInit[0], fZPlane, 0., dpy);
1938 fSeg[0]->MorePads(); fSeg[0]->NextPad())
1939 {
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];
1944 icount++;
1945 }
1946
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]);
1949
1950
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);
1957
1958
1959// Find save upper and lower limits
1960
1961 icount=0;
1962
1963 for (fSeg[1]->FirstPad(fXInit[1], fYInit[1], fZPlane, dpx, 0);
1964 fSeg[1]->MorePads(); fSeg[1]->NextPad())
1965 {
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];
1970 icount++;
1971 }
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]);
1974
1975 icount=0;
1976
1977 for (fSeg[0]->FirstPad(fXInit[1], fYInit[1], fZPlane, 0, dpy);
1978 fSeg[0]-> MorePads(); fSeg[0]->NextPad())
1979 {
1980 ix=fSeg[0]->Ix(); iy=fSeg[0]->Iy();
1981// if (fHitMap[0]->TestHit(ix, iy) != kEmpty) continue;
1982
1983 fSeg[0]->GetPadC(ix,iy,xdum,upper[3],zdum);
1984 if (icount ==0) lower[3]=upper[3];
1985 icount++;
1986
1987 }
1988 if (lower[3]>upper[3]) {xdum=lower[3]; lower[3]=upper[3]; upper[3]=xdum;}
1989
1990// vstart[3] = 0.5*(lower[3]+upper[3]);
1991
1992 lower[4]=0.;
1993 upper[4]=1.;
1994 lower[5]=0.;
1995 upper[5]=1.;
1996
1997// step sizes
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);
2007 if (fDebugLevel)
2008 clusterInput.Fitter()->mnexcm("SET NOW", arglist, 0, ierflag);
2009 clusterInput.Fitter()->mnexcm("SET OUT", arglist, 0, ierflag);
2010 arglist[0]= -1;
2011 arglist[1]= 0;
2012
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
2017 TString chname;
2018 Double_t epxz, b1, b2;
2019 Int_t ierflg;
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);
2026
2027 Double_t fmin, fedm, errdef;
2028 Int_t npari, nparx, istat;
2029
2030 clusterInput.Fitter()->mnstat(fmin, fedm, errdef, npari, nparx, istat);
2031 fFitStat=istat;
2032
2033 fChi2[0]=fmin;
2034 fChi2[1]=fmin;
2035 return fmin;
2036}
2037
2038void AliMUONClusterFinderVS::Split(AliMUONRawCluster* c)
2039{
2040//
2041// One cluster for each maximum
2042//
2043 Int_t i, j, cath;
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];
2051
2052 if (fNPeaks == 0) {
2053 cnew.fNcluster[0]=-1;
2054 cnew.fNcluster[1]=fNRawClusters;
2055 } else {
2056 cnew.fNcluster[0]=fNPeaks;
2057 cnew.fNcluster[1]=0;
2058 }
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;
2063 if (j==0) {
2064 cnew.fQ[cath]=Int_t(clusterInput.TotalCharge(cath)*fQrFit[cath]);
2065 } else {
2066 cnew.fQ[cath]=Int_t(clusterInput.TotalCharge(cath)*(1-fQrFit[cath]));
2067 }
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]++;
2077 }
2078 FillCluster(&cnew,0,cath);
2079 } // cathode loop
2080
2081 cnew.fClusterType=cnew.PhysicsContribution();
2082 if (cnew.fQ[0]>0 && cnew.fQ[1]>0) AddRawCluster(cnew);
2083 fNPeaks++;
2084 }
2085}
2086
2087
2088//
2089// Minimisation functions
2090// Single Mathieson
2091void fcnS1(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag)
2092{
2093 AliMUONClusterInput& clusterInput = *(AliMUONClusterInput::Instance());
2094 Int_t i;
2095 Float_t delta;
2096 Float_t chisq=0;
2097 Float_t qcont=0;
2098 Float_t qtot=0;
2099
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);
2103 delta=(q0-q1)/q0;
2104 chisq+=delta*delta;
2105 qcont+=q1;
2106 qtot+=q0;
2107 }
2108 f=chisq;
2109}
2110
2111void fcnCombiS1(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag)
2112{
2113 AliMUONClusterInput& clusterInput = *(AliMUONClusterInput::Instance());
2114 Int_t i, cath;
2115 Float_t delta;
2116 Float_t chisq=0;
2117 Float_t qcont=0;
2118 Float_t qtot=0;
2119
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);
2124 delta=(q0-q1)/q0;
2125 chisq+=delta*delta;
2126 qcont+=q1;
2127 qtot+=q0;
2128 }
2129 }
2130 f=chisq;
2131}
2132
2133// Double Mathieson
2134void fcnS2(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag)
2135{
2136 AliMUONClusterInput& clusterInput = *(AliMUONClusterInput::Instance());
2137 Int_t i;
2138 Float_t delta;
2139 Float_t chisq=0;
2140 Float_t qcont=0;
2141 Float_t qtot=0;
2142
2143 for (i=0; i<clusterInput.Nmul(0); i++) {
2144
2145 Float_t q0=clusterInput.Charge(i,0);
2146 Float_t q1=clusterInput.DiscrChargeS2(i,par);
2147 delta=(q0-q1)/q0;
2148 chisq+=delta*delta;
2149 qcont+=q1;
2150 qtot+=q0;
2151 }
2152 f=chisq;
2153}
2154
2155// Double Mathieson
2156void fcnCombiS2(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag)
2157{
2158 AliMUONClusterInput& clusterInput = *(AliMUONClusterInput::Instance());
2159 Int_t i, cath;
2160 Float_t delta;
2161 Float_t chisq=0;
2162 Float_t qcont=0;
2163 Float_t qtot=0;
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);
2168 delta=(q0-q1)/q0;
2169 chisq+=delta*delta;
2170 qcont+=q1;
2171 qtot+=q0;
2172 }
2173 }
2174 f=chisq;
2175}
2176
2177void AliMUONClusterFinderVS::AddRawCluster(const AliMUONRawCluster c)
2178{
2179 //
2180 // Add a raw cluster copy to the list
2181 //
2182 AliMUON *pMUON=(AliMUON*)gAlice->GetModule("MUON");
2183 pMUON->AddRawCluster(fInput->Chamber(),c);
2184 fNRawClusters++;
2185// if (fDebugLevel)
2186 fprintf(stderr,"\nfNRawClusters %d\n",fNRawClusters);
2187}
2188
2189Bool_t AliMUONClusterFinderVS::TestTrack(Int_t t) {
2190// Test if track was user selected
2191 if (fTrack[0]==-1 || fTrack[1]==-1) {
2192 return kTRUE;
2193 } else if (t==fTrack[0] || t==fTrack[1]) {
2194 return kTRUE;
2195 } else {
2196 return kFALSE;
2197 }
2198}
2199
2200AliMUONClusterFinderVS& AliMUONClusterFinderVS
2201::operator = (const AliMUONClusterFinderVS& rhs)
2202{
2203// Dummy assignment operator
2204 return *this;
2205}
2206
2207
2208
2209
2210
2211
2212
2213
2214