]> 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
AliMUONSegmentation and AliMUONHitMap have been replaced by AliSegmentation and AliHi...
[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.7 2000/06/28 15:16:35 morsch
18(1) Client code adapted to new method signatures in AliMUONSegmentation (see comments there)
19to allow development of slat-muon chamber simulation and reconstruction code in the MUON
20framework. The changes should have no side effects (mostly dummy arguments).
21(2) Hit disintegration uses 3-dim hit coordinates to allow simulation
22of chambers with overlapping modules (MakePadHits, Disintegration).
23
24Revision 1.6 2000/06/28 12:19:18 morsch
25More consequent seperation of global input data services (AliMUONClusterInput singleton) and the
26cluster and hit reconstruction algorithms in AliMUONClusterFinderVS.
27AliMUONClusterFinderVS becomes the base class for clustering and hit reconstruction.
28It requires two cathode planes. Small modifications in the code will make it usable for
29one cathode plane and, hence, more general (for test beam data).
30AliMUONClusterFinder is now obsolete.
31
32Revision 1.5 2000/06/28 08:06:10 morsch
33Avoid global variables in AliMUONClusterFinderVS by seperating the input data for the fit from the
34algorithmic part of the class. Input data resides inside the AliMUONClusterInput singleton.
35It also naturally takes care of the TMinuit instance.
36
37Revision 1.4 2000/06/27 16:18:47 gosset
38Finally correct implementation of xm, ym, ixm, iym sizes
39when at least three local maxima on cathode 1 or on cathode 2
40
41Revision 1.3 2000/06/22 14:02:45 morsch
42Parameterised size of xm[], ym[], ixm[], iym[] correctly implemented (PH)
43Some HP scope problems corrected (PH)
44
45Revision 1.2 2000/06/15 07:58:48 morsch
46Code from MUON-dev joined
47
48Revision 1.1.2.3 2000/06/09 21:58:33 morsch
49Most coding rule violations corrected.
50
51Revision 1.1.2.2 2000/02/15 08:33:52 morsch
52Error in calculation of contribution map for double clusters (Split method) corrected (A.M.)
53Error in determination of track list for double cluster (FillCluster method) corrected (A.M.)
54Revised and extended SplitByLocalMaxima method (Isabelle Chevrot):
55 - For clusters with more than 2 maxima on one of the cathode planes all valid
56 combinations of maxima on the two cathodes are preserved. The position of the maxima is
57 taken as the hit position.
58 - New FillCluster method with 2 arguments to find tracks associated to the clusters
59 defined above added. (Method destinction by argument list not very elegant in this case,
60 should be revides (A.M.)
61 - Bug in if-statement to handle maximum 1 maximum per plane corrected
62 - Two cluster per cathode but only 1 combination valid is handled.
63 - More rigerous treatment of 1-2 and 2-1 combinations of maxima.
64
65*/
66
67#include "AliMUONClusterFinderVS.h"
68#include "AliMUONDigit.h"
69#include "AliMUONRawCluster.h"
70#include "AliSegmentation.h"
71#include "AliMUONResponse.h"
72#include "AliMUONHitMapA1.h"
73#include "AliRun.h"
74#include "AliMUON.h"
75
76#include <TTree.h>
77#include <TCanvas.h>
78#include <TH1.h>
79#include <TPad.h>
80#include <TGraph.h>
81#include <TPostScript.h>
82#include <TMinuit.h>
83#include <stdio.h>
84#include <iostream.h>
85
86//_____________________________________________________________________
87// This function is minimized in the double-Mathieson fit
88void fcnS2(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag);
89void fcnS1(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag);
90void fcnCombiS1(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag);
91void fcnCombiS2(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag);
92
93ClassImp(AliMUONClusterFinderVS)
94
95 AliMUONClusterFinderVS::AliMUONClusterFinderVS()
96{
97// Default constructor
98 fInput=AliMUONClusterInput::Instance();
99 fHitMap[0] = 0;
100 fHitMap[1] = 0;
101 fTrack[0]=fTrack[1]=-1;
102}
103
104AliMUONClusterFinderVS::AliMUONClusterFinderVS(
105 const AliMUONClusterFinderVS & clusterFinder)
106{
107// Dummy copy Constructor
108 ;
109}
110
111void AliMUONClusterFinderVS::Decluster(AliMUONRawCluster *cluster)
112{
113// Decluster by local maxima
114 SplitByLocalMaxima(cluster);
115}
116
117void AliMUONClusterFinderVS::SplitByLocalMaxima(AliMUONRawCluster *c)
118{
119// Split complex cluster by local maxima
120
121 Int_t cath, i;
122
123 fInput->SetCluster(c);
124
125 fMul[0]=c->fMultiplicity[0];
126 fMul[1]=c->fMultiplicity[1];
127
128//
129// dump digit information into arrays
130//
131
132 Float_t qtot, zdum;
133
134 for (cath=0; cath<2; cath++) {
135 qtot=0;
136 for (i=0; i<fMul[cath]; i++)
137 {
138 // pointer to digit
139 fDig[i][cath]=fInput->Digit(cath, c->fIndexMap[i][cath]);
140 // pad coordinates
141 fIx[i][cath]= fDig[i][cath]->fPadX;
142 fIy[i][cath]= fDig[i][cath]->fPadY;
143 // pad charge
144 fQ[i][cath] = fDig[i][cath]->fSignal;
145 // pad centre coordinates
146 fInput->Segmentation(cath)->
147 GetPadC(fIx[i][cath], fIy[i][cath], fX[i][cath], fY[i][cath], zdum);
148 } // loop over cluster digits
149 } // loop over cathodes
150
151
152 FindLocalMaxima(c);
153
154//
155// Initialise and perform mathieson fits
156 Float_t chi2, oldchi2;
157// ++++++++++++++++++*************+++++++++++++++++++++
158// (1) No more than one local maximum per cathode plane
159// +++++++++++++++++++++++++++++++*************++++++++
160 if ((fNLocal[0]==1 && (fNLocal[1]==0 || fNLocal[1]==1)) ||
161 (fNLocal[0]==0 && fNLocal[1]==1)) {
162
163// Perform combined single Mathieson fit
164// Initial values for coordinates (x,y)
165
166 // One local maximum on cathodes 1 and 2 (X->cathode 2, Y->cathode 1)
167 if (fNLocal[0]==1 && fNLocal[1]==1) {
168 fXInit[0]=c->fX[1];
169 fYInit[0]=c->fY[0];
170 // One local maximum on cathode 1 (X,Y->cathode 1)
171 } else if (fNLocal[0]==1) {
172 fXInit[0]=c->fX[0];
173 fYInit[0]=c->fY[0];
174 // One local maximum on cathode 2 (X,Y->cathode 2)
175 } else {
176 fXInit[0]=c->fX[1];
177 fYInit[0]=c->fY[1];
178 }
179 fprintf(stderr,"\n cas (1) CombiSingleMathiesonFit(c)\n");
180 chi2=CombiSingleMathiesonFit(c);
181// Int_t ndf = fgNbins[0]+fgNbins[1]-2;
182// Float_t prob = TMath::Prob(Double_t(chi2),ndf);
183// prob1->Fill(prob);
184// chi2_1->Fill(chi2);
185 oldchi2=chi2;
186 fprintf(stderr," chi2 %f ",chi2);
187
188 c->fX[0]=fXFit[0];
189 c->fY[0]=fYFit[0];
190
191 c->fX[1]=fXFit[0];
192 c->fY[1]=fYFit[0];
193 c->fChi2[0]=chi2;
194 c->fChi2[1]=chi2;
195 c->fX[0]=fInput->Segmentation(0)->GetAnod(c->fX[0]);
196 c->fX[1]=fInput->Segmentation(1)->GetAnod(c->fX[1]);
197
198// If reasonable chi^2 add result to the list of rawclusters
199 // if (chi2 < 50) {
200 if (chi2 < 0.3) {
201 AddRawCluster(*c);
202// If not try combined double Mathieson Fit
203 } else {
204 fprintf(stderr," MAUVAIS CHI2 !!!\n");
205 if (fNLocal[0]==1 && fNLocal[1]==1) {
206 fXInit[0]=fX[fIndLocal[0][1]][1];
207 fYInit[0]=fY[fIndLocal[0][0]][0];
208 fXInit[1]=fX[fIndLocal[0][1]][1];
209 fYInit[1]=fY[fIndLocal[0][0]][0];
210 } else if (fNLocal[0]==1) {
211 fXInit[0]=fX[fIndLocal[0][0]][0];
212 fYInit[0]=fY[fIndLocal[0][0]][0];
213 fXInit[1]=fX[fIndLocal[0][0]][0];
214 fYInit[1]=fY[fIndLocal[0][0]][0];
215 } else {
216 fXInit[0]=fX[fIndLocal[0][1]][1];
217 fYInit[0]=fY[fIndLocal[0][1]][1];
218 fXInit[1]=fX[fIndLocal[0][1]][1];
219 fYInit[1]=fY[fIndLocal[0][1]][1];
220 }
221
222// Initial value for charge ratios
223 fQrInit[0]=0.5;
224 fQrInit[1]=0.5;
225 fprintf(stderr,"\n cas (1) CombiDoubleMathiesonFit(c)\n");
226 chi2=CombiDoubleMathiesonFit(c);
227// Int_t ndf = fgNbins[0]+fgNbins[1]-6;
228// Float_t prob = TMath::Prob(chi2,ndf);
229// prob2->Fill(prob);
230// chi2_2->Fill(chi2);
231
232// Was this any better ??
233 fprintf(stderr," Old and new chi2 %f %f ", oldchi2, chi2);
234 if (fFitStat!=0 && chi2>0 && (2.*chi2 < oldchi2)) {
235 fprintf(stderr," Split\n");
236 // Split cluster into two according to fit result
237 Split(c);
238 } else {
239 fprintf(stderr," Don't Split\n");
240 // Don't split
241 AddRawCluster(*c);
242 }
243 }
244
245// +++++++++++++++++++++++++++++++++++++++
246// (2) Two local maxima per cathode plane
247// +++++++++++++++++++++++++++++++++++++++
248 } else if (fNLocal[0]==2 && fNLocal[1]==2) {
249//
250// Let's look for ghosts first
251//
252 Float_t xm[4][2], ym[4][2];
253 Float_t dpx, dpy, dx, dy;
254 Int_t ixm[4][2], iym[4][2];
255 Int_t isec, im1, im2, ico;
256//
257// Form the 2x2 combinations
258// 0-0, 0-1, 1-0, 1-1
259 ico=0;
260 for (im1=0; im1<2; im1++) {
261 for (im2=0; im2<2; im2++) {
262 xm[ico][0]=fX[fIndLocal[im1][0]][0];
263 ym[ico][0]=fY[fIndLocal[im1][0]][0];
264 xm[ico][1]=fX[fIndLocal[im2][1]][1];
265 ym[ico][1]=fY[fIndLocal[im2][1]][1];
266
267 ixm[ico][0]=fIx[fIndLocal[im1][0]][0];
268 iym[ico][0]=fIy[fIndLocal[im1][0]][0];
269 ixm[ico][1]=fIx[fIndLocal[im2][1]][1];
270 iym[ico][1]=fIy[fIndLocal[im2][1]][1];
271 ico++;
272 }
273 }
274// ico = 0 : first local maximum on cathodes 1 and 2
275// ico = 1 : fisrt local maximum on cathode 1 and second on cathode 2
276// ico = 2 : second local maximum on cathode 1 and first on cathode 1
277// ico = 3 : second local maximum on cathodes 1 and 2
278
279// Analyse the combinations and keep those that are possible !
280// For each combination check consistency in x and y
281 Int_t iacc;
282 Bool_t accepted[4];
283 iacc=0;
284
285 for (ico=0; ico<4; ico++) {
286 accepted[ico]=kFALSE;
287// cathode one: x-coordinate
288 isec=fInput->Segmentation(0)->Sector(ixm[ico][0], iym[ico][0]);
289 dpx=fInput->Segmentation(0)->Dpx(isec)/2.;
290 dx=TMath::Abs(xm[ico][0]-xm[ico][1]);
291// cathode two: y-coordinate
292 isec=fInput->Segmentation(1)->Sector(ixm[ico][1], iym[ico][1]);
293 dpy=fInput->Segmentation(1)->Dpy(isec)/2.;
294 dy=TMath::Abs(ym[ico][0]-ym[ico][1]);
295// printf("\n %i %f %f %f %f \n", ico, ym[ico][0], ym[ico][1], dy, dpy );
296 if ((dx <= dpx) && (dy <= dpy)) {
297 // consistent
298 accepted[ico]=kTRUE;
299 iacc++;
300 } else {
301 // reject
302 accepted[ico]=kFALSE;
303 }
304 }
305
306 if (iacc==2) {
307 fprintf(stderr,"\n iacc=2: No problem ! \n");
308 } else if (iacc==4) {
309 fprintf(stderr,"\n iacc=4: Ok, but ghost problem !!! \n");
310 } else if (iacc==0) {
311 fprintf(stderr,"\n iacc=0: I don't know what to do with this !!!!!!!!! \n");
312 }
313
314// Initial value for charge ratios
315 fQrInit[0]=Float_t(fQ[fIndLocal[0][0]][0])/
316 Float_t(fQ[fIndLocal[0][0]][0]+fQ[fIndLocal[1][0]][0]);
317 fQrInit[1]=Float_t(fQ[fIndLocal[0][1]][1])/
318 Float_t(fQ[fIndLocal[0][1]][1]+fQ[fIndLocal[1][1]][1]);
319
320// ******* iacc = 0 *******
321// No combinations found between the 2 cathodes
322// We keep the center of gravity of the cluster
323 if (iacc==0) {
324 AddRawCluster(*c);
325 }
326
327// ******* iacc = 1 *******
328// Only one combination found between the 2 cathodes
329 if (iacc==1) {
330
331// Initial values for the 2 maxima (x,y)
332
333// 1 maximum is initialised with the maximum of the combination found (X->cathode 2, Y->cathode 1)
334// 1 maximum is initialised with the other maximum of the first cathode
335 if (accepted[0]){
336 fprintf(stderr,"ico=0\n");
337 fXInit[0]=xm[0][1];
338 fYInit[0]=ym[0][0];
339 fXInit[1]=xm[3][0];
340 fYInit[1]=ym[3][0];
341 } else if (accepted[1]){
342 fprintf(stderr,"ico=1\n");
343 fXInit[0]=xm[1][1];
344 fYInit[0]=ym[1][0];
345 fXInit[1]=xm[2][0];
346 fYInit[1]=ym[2][0];
347 } else if (accepted[2]){
348 fprintf(stderr,"ico=2\n");
349 fXInit[0]=xm[2][1];
350 fYInit[0]=ym[2][0];
351 fXInit[1]=xm[1][0];
352 fYInit[1]=ym[1][0];
353 } else if (accepted[3]){
354 fprintf(stderr,"ico=3\n");
355 fXInit[0]=xm[3][1];
356 fYInit[0]=ym[3][0];
357 fXInit[1]=xm[0][0];
358 fYInit[1]=ym[0][0];
359 }
360 fprintf(stderr,"\n cas (2) CombiDoubleMathiesonFit(c)\n");
361 chi2=CombiDoubleMathiesonFit(c);
362// Int_t ndf = fgNbins[0]+fgNbins[1]-6;
363// Float_t prob = TMath::Prob(chi2,ndf);
364// prob2->Fill(prob);
365// chi2_2->Fill(chi2);
366 fprintf(stderr," chi2 %f\n",chi2);
367
368// If reasonable chi^2 add result to the list of rawclusters
369 if (chi2<10) {
370 Split(c);
371
372 } else {
373// 1 maximum is initialised with the maximum of the combination found (X->cathode 2, Y->cathode 1)
374// 1 maximum is initialised with the other maximum of the second cathode
375 if (accepted[0]){
376 fprintf(stderr,"ico=0\n");
377 fXInit[0]=xm[0][1];
378 fYInit[0]=ym[0][0];
379 fXInit[1]=xm[3][1];
380 fYInit[1]=ym[3][1];
381 } else if (accepted[1]){
382 fprintf(stderr,"ico=1\n");
383 fXInit[0]=xm[1][1];
384 fYInit[0]=ym[1][0];
385 fXInit[1]=xm[2][1];
386 fYInit[1]=ym[2][1];
387 } else if (accepted[2]){
388 fprintf(stderr,"ico=2\n");
389 fXInit[0]=xm[2][1];
390 fYInit[0]=ym[2][0];
391 fXInit[1]=xm[1][1];
392 fYInit[1]=ym[1][1];
393 } else if (accepted[3]){
394 fprintf(stderr,"ico=3\n");
395 fXInit[0]=xm[3][1];
396 fYInit[0]=ym[3][0];
397 fXInit[1]=xm[0][1];
398 fYInit[1]=ym[0][1];
399 }
400 fprintf(stderr,"\n cas (2) CombiDoubleMathiesonFit(c)\n");
401 chi2=CombiDoubleMathiesonFit(c);
402// Int_t ndf = fgNbins[0]+fgNbins[1]-6;
403// Float_t prob = TMath::Prob(chi2,ndf);
404// prob2->Fill(prob);
405// chi2_2->Fill(chi2);
406 fprintf(stderr," chi2 %f\n",chi2);
407
408// If reasonable chi^2 add result to the list of rawclusters
409 if (chi2<10) {
410 Split(c);
411 } else {
412//We keep only the combination found (X->cathode 2, Y->cathode 1)
413 for (Int_t ico=0; ico<2; ico++) {
414 if (accepted[ico]) {
415 AliMUONRawCluster cnew;
416 Int_t cath;
417 for (cath=0; cath<2; cath++) {
418 cnew.fX[cath]=Float_t(xm[ico][1]);
419 cnew.fY[cath]=Float_t(ym[ico][0]);
420 cnew.fMultiplicity[cath]=c->fMultiplicity[cath];
421 for (i=0; i<fMul[cath]; i++) {
422 cnew.fIndexMap[i][cath]=c->fIndexMap[i][cath];
423 fInput->Segmentation(cath)->SetPad(fIx[i][cath], fIy[i][cath]);
424 }
425 fprintf(stderr,"\nRawCluster %d cath %d\n",ico,cath);
426 fprintf(stderr,"mult_av %d\n",c->fMultiplicity[cath]);
427 FillCluster(&cnew,cath);
428 }
429 cnew.fClusterType=cnew.PhysicsContribution();
430 AddRawCluster(cnew);
431 fNPeaks++;
432 }
433 }
434 }
435 }
436 }
437
438// ******* iacc = 2 *******
439// Two combinations found between the 2 cathodes
440 if (iacc==2) {
441
442// Was the same maximum taken twice
443 if ((accepted[0]&&accepted[1]) || (accepted[2]&&accepted[3])) {
444 fprintf(stderr,"\n Maximum taken twice !!!\n");
445
446// Have a try !! with that
447 if (accepted[0]&&accepted[3]) {
448 fXInit[0]=xm[0][1];
449 fYInit[0]=ym[0][0];
450 fXInit[1]=xm[1][1];
451 fYInit[1]=ym[1][0];
452 } else {
453 fXInit[0]=xm[2][1];
454 fYInit[0]=ym[2][0];
455 fXInit[1]=xm[3][1];
456 fYInit[1]=ym[3][0];
457 }
458 fprintf(stderr,"\n cas (2) CombiDoubleMathiesonFit(c)\n");
459 chi2=CombiDoubleMathiesonFit(c);
460// Int_t ndf = fgNbins[0]+fgNbins[1]-6;
461// Float_t prob = TMath::Prob(chi2,ndf);
462// prob2->Fill(prob);
463// chi2_2->Fill(chi2);
464 Split(c);
465
466 } else {
467// No ghosts ! No Problems ! - Perform one fit only !
468 if (accepted[0]&&accepted[3]) {
469 fXInit[0]=xm[0][1];
470 fYInit[0]=ym[0][0];
471 fXInit[1]=xm[3][1];
472 fYInit[1]=ym[3][0];
473 } else {
474 fXInit[0]=xm[1][1];
475 fYInit[0]=ym[1][0];
476 fXInit[1]=xm[2][1];
477 fYInit[1]=ym[2][0];
478 }
479 fprintf(stderr,"\n cas (2) CombiDoubleMathiesonFit(c)\n");
480 chi2=CombiDoubleMathiesonFit(c);
481// Int_t ndf = fgNbins[0]+fgNbins[1]-6;
482// Float_t prob = TMath::Prob(chi2,ndf);
483// prob2->Fill(prob);
484// chi2_2->Fill(chi2);
485 fprintf(stderr," chi2 %f\n",chi2);
486 Split(c);
487 }
488
489// ******* iacc = 4 *******
490// Four combinations found between the 2 cathodes
491// Ghost !!
492 } else if (iacc==4) {
493// Perform fits for the two possibilities !!
494 fXInit[0]=xm[0][1];
495 fYInit[0]=ym[0][0];
496 fXInit[1]=xm[3][1];
497 fYInit[1]=ym[3][0];
498 fprintf(stderr,"\n cas (2) CombiDoubleMathiesonFit(c)\n");
499 chi2=CombiDoubleMathiesonFit(c);
500// Int_t ndf = fgNbins[0]+fgNbins[1]-6;
501// Float_t prob = TMath::Prob(chi2,ndf);
502// prob2->Fill(prob);
503// chi2_2->Fill(chi2);
504 fprintf(stderr," chi2 %f\n",chi2);
505 Split(c);
506 fXInit[0]=xm[1][1];
507 fYInit[0]=ym[1][0];
508 fXInit[1]=xm[2][1];
509 fYInit[1]=ym[2][0];
510 fprintf(stderr,"\n cas (2) CombiDoubleMathiesonFit(c)\n");
511 chi2=CombiDoubleMathiesonFit(c);
512// ndf = fgNbins[0]+fgNbins[1]-6;
513// prob = TMath::Prob(chi2,ndf);
514// prob2->Fill(prob);
515// chi2_2->Fill(chi2);
516 fprintf(stderr," chi2 %f\n",chi2);
517 Split(c);
518 }
519
520 } else if (fNLocal[0]==2 && fNLocal[1]==1) {
521// ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
522// (3) Two local maxima on cathode 1 and one maximum on cathode 2
523// ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
524//
525 Float_t xm[4][2], ym[4][2];
526 Float_t dpx, dpy, dx, dy;
527 Int_t ixm[4][2], iym[4][2];
528 Int_t isec, im1, ico;
529//
530// Form the 2x2 combinations
531// 0-0, 0-1, 1-0, 1-1
532 ico=0;
533 for (im1=0; im1<2; im1++) {
534 xm[ico][0]=fX[fIndLocal[im1][0]][0];
535 ym[ico][0]=fY[fIndLocal[im1][0]][0];
536 xm[ico][1]=fX[fIndLocal[0][1]][1];
537 ym[ico][1]=fY[fIndLocal[0][1]][1];
538
539 ixm[ico][0]=fIx[fIndLocal[im1][0]][0];
540 iym[ico][0]=fIy[fIndLocal[im1][0]][0];
541 ixm[ico][1]=fIx[fIndLocal[0][1]][1];
542 iym[ico][1]=fIy[fIndLocal[0][1]][1];
543 ico++;
544 }
545// ico = 0 : first local maximum on cathodes 1 and 2
546// ico = 1 : second local maximum on cathode 1 and first on cathode 2
547
548// Analyse the combinations and keep those that are possible !
549// For each combination check consistency in x and y
550 Int_t iacc;
551 Bool_t accepted[4];
552 iacc=0;
553
554 for (ico=0; ico<2; ico++) {
555 accepted[ico]=kFALSE;
556 isec=fInput->Segmentation(0)->Sector(ixm[ico][0], iym[ico][0]);
557 dpx=fInput->Segmentation(0)->Dpx(isec)/2.;
558 dx=TMath::Abs(xm[ico][0]-xm[ico][1]);
559 isec=fInput->Segmentation(1)->Sector(ixm[ico][1], iym[ico][1]);
560 dpy=fInput->Segmentation(1)->Dpy(isec)/2.;
561 dy=TMath::Abs(ym[ico][0]-ym[ico][1]);
562// printf("\n %i %f %f %f %f \n", ico, ym[ico][0], ym[ico][1], dy, dpy );
563 if ((dx <= dpx) && (dy <= dpy)) {
564 // consistent
565 accepted[ico]=kTRUE;
566 iacc++;
567 } else {
568 // reject
569 accepted[ico]=kFALSE;
570 }
571 }
572
573 Float_t chi21 = 100;
574 Float_t chi22 = 100;
575
576 if (accepted[0]) {
577 fXInit[0]=xm[0][1];
578 fYInit[0]=ym[0][0];
579 fXInit[1]=xm[1][0];
580 fYInit[1]=ym[1][0];
581 chi21=CombiDoubleMathiesonFit(c);
582// Int_t ndf = fgNbins[0]+fgNbins[1]-6;
583// Float_t prob = TMath::Prob(chi2,ndf);
584// prob2->Fill(prob);
585// chi2_2->Fill(chi21);
586 fprintf(stderr," chi2 %f\n",chi21);
587 if (chi21<10) Split(c);
588 } else if (accepted[1]) {
589 fXInit[0]=xm[1][1];
590 fYInit[0]=ym[1][0];
591 fXInit[1]=xm[0][0];
592 fYInit[1]=ym[0][0];
593 chi22=CombiDoubleMathiesonFit(c);
594// Int_t ndf = fgNbins[0]+fgNbins[1]-6;
595// Float_t prob = TMath::Prob(chi2,ndf);
596// prob2->Fill(prob);
597// chi2_2->Fill(chi22);
598 fprintf(stderr," chi2 %f\n",chi22);
599 if (chi22<10) Split(c);
600 }
601
602 if (chi21 > 10 && chi22 > 10) {
603// We keep only the combination found (X->cathode 2, Y->cathode 1)
604 for (Int_t ico=0; ico<2; ico++) {
605 if (accepted[ico]) {
606 AliMUONRawCluster cnew;
607 Int_t cath;
608 for (cath=0; cath<2; cath++) {
609 cnew.fX[cath]=Float_t(xm[ico][1]);
610 cnew.fY[cath]=Float_t(ym[ico][0]);
611 cnew.fMultiplicity[cath]=c->fMultiplicity[cath];
612 for (i=0; i<fMul[cath]; i++) {
613 cnew.fIndexMap[i][cath]=c->fIndexMap[i][cath];
614 fInput->Segmentation(cath)->SetPad(fIx[i][cath], fIy[i][cath]);
615 }
616 fprintf(stderr,"\nRawCluster %d cath %d\n",ico,cath);
617 fprintf(stderr,"mult_av %d\n",c->fMultiplicity[cath]);
618 FillCluster(&cnew,cath);
619 }
620 cnew.fClusterType=cnew.PhysicsContribution();
621 AddRawCluster(cnew);
622 fNPeaks++;
623 }
624 }
625 }
626
627// +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
628// (3') One local maximum on cathode 1 and two maxima on cathode 2
629// +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
630 } else if (fNLocal[0]==1 && fNLocal[1]==2) {
631
632 Float_t xm[4][2], ym[4][2];
633 Float_t dpx, dpy, dx, dy;
634 Int_t ixm[4][2], iym[4][2];
635 Int_t isec, im1, ico;
636//
637// Form the 2x2 combinations
638// 0-0, 0-1, 1-0, 1-1
639 ico=0;
640 for (im1=0; im1<2; im1++) {
641 xm[ico][0]=fX[fIndLocal[0][0]][0];
642 ym[ico][0]=fY[fIndLocal[0][0]][0];
643 xm[ico][1]=fX[fIndLocal[im1][1]][1];
644 ym[ico][1]=fY[fIndLocal[im1][1]][1];
645
646 ixm[ico][0]=fIx[fIndLocal[0][0]][0];
647 iym[ico][0]=fIy[fIndLocal[0][0]][0];
648 ixm[ico][1]=fIx[fIndLocal[im1][1]][1];
649 iym[ico][1]=fIy[fIndLocal[im1][1]][1];
650 ico++;
651 }
652// ico = 0 : first local maximum on cathodes 1 and 2
653// ico = 1 : first local maximum on cathode 1 and second on cathode 2
654
655// Analyse the combinations and keep those that are possible !
656// For each combination check consistency in x and y
657 Int_t iacc;
658 Bool_t accepted[4];
659 iacc=0;
660
661 for (ico=0; ico<2; ico++) {
662 accepted[ico]=kFALSE;
663 isec=fInput->Segmentation(0)->Sector(ixm[ico][0], iym[ico][0]);
664 dpx=fInput->Segmentation(0)->Dpx(isec)/2.;
665 dx=TMath::Abs(xm[ico][0]-xm[ico][1]);
666 isec=fInput->Segmentation(1)->Sector(ixm[ico][1], iym[ico][1]);
667 dpy=fInput->Segmentation(1)->Dpy(isec)/2.;
668 dy=TMath::Abs(ym[ico][0]-ym[ico][1]);
669// printf("\n %i %f %f %f %f \n", ico, ym[ico][0], ym[ico][1], dy, dpy );
670 if ((dx <= dpx) && (dy <= dpy)) {
671 // consistent
672 accepted[ico]=kTRUE;
673 fprintf(stderr,"ico %d\n",ico);
674 iacc++;
675 } else {
676 // reject
677 accepted[ico]=kFALSE;
678 }
679 }
680
681 Float_t chi21 = 100;
682 Float_t chi22 = 100;
683
684 if (accepted[0]) {
685 fXInit[0]=xm[0][0];
686 fYInit[0]=ym[0][1];
687 fXInit[1]=xm[1][1];
688 fYInit[1]=ym[1][1];
689 chi21=CombiDoubleMathiesonFit(c);
690// Int_t ndf = fgNbins[0]+fgNbins[1]-6;
691// Float_t prob = TMath::Prob(chi2,ndf);
692// prob2->Fill(prob);
693// chi2_2->Fill(chi21);
694 fprintf(stderr," chi2 %f\n",chi21);
695 if (chi21<10) Split(c);
696 } else if (accepted[1]) {
697 fXInit[0]=xm[1][0];
698 fYInit[0]=ym[1][1];
699 fXInit[1]=xm[0][1];
700 fYInit[1]=ym[0][1];
701 chi22=CombiDoubleMathiesonFit(c);
702// Int_t ndf = fgNbins[0]+fgNbins[1]-6;
703// Float_t prob = TMath::Prob(chi2,ndf);
704// prob2->Fill(prob);
705// chi2_2->Fill(chi22);
706 fprintf(stderr," chi2 %f\n",chi22);
707 if (chi22<10) Split(c);
708 }
709
710 if (chi21 > 10 && chi22 > 10) {
711//We keep only the combination found (X->cathode 2, Y->cathode 1)
712 for (Int_t ico=0; ico<2; ico++) {
713 if (accepted[ico]) {
714 AliMUONRawCluster cnew;
715 Int_t cath;
716 for (cath=0; cath<2; cath++) {
717 cnew.fX[cath]=Float_t(xm[ico][1]);
718 cnew.fY[cath]=Float_t(ym[ico][0]);
719 cnew.fMultiplicity[cath]=c->fMultiplicity[cath];
720 for (i=0; i<fMul[cath]; i++) {
721 cnew.fIndexMap[i][cath]=c->fIndexMap[i][cath];
722 fInput->Segmentation(cath)->SetPad(fIx[i][cath], fIy[i][cath]);
723 }
724 fprintf(stderr,"\nRawCluster %d cath %d\n",ico,cath);
725 fprintf(stderr,"mult_av %d\n",c->fMultiplicity[cath]);
726 FillCluster(&cnew,cath);
727 }
728 cnew.fClusterType=cnew.PhysicsContribution();
729 AddRawCluster(cnew);
730 fNPeaks++;
731 }
732 }
733 }
734
735// ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
736// (4) At least three local maxima on cathode 1 or on cathode 2
737// ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
738 } else if (fNLocal[0]>2 || fNLocal[1]>2) {
739
740 Int_t param = fNLocal[0]*fNLocal[1];
741 Int_t ii;
742
743 Float_t ** xm = new Float_t * [param];
744 for (ii=0; ii<param; ii++) xm[ii]=new Float_t [2];
745 Float_t ** ym = new Float_t * [param];
746 for (ii=0; ii<param; ii++) ym[ii]=new Float_t [2];
747 Int_t ** ixm = new Int_t * [param];
748 for (ii=0; ii<param; ii++) ixm[ii]=new Int_t [2];
749 Int_t ** iym = new Int_t * [param];
750 for (ii=0; ii<param; ii++) iym[ii]=new Int_t [2];
751
752 Int_t isec, ico;
753 Float_t dpx, dpy, dx, dy;
754
755 ico=0;
756 for (Int_t im1=0; im1<fNLocal[0]; im1++) {
757 for (Int_t im2=0; im2<fNLocal[1]; im2++) {
758 xm[ico][0]=fX[fIndLocal[im1][0]][0];
759 ym[ico][0]=fY[fIndLocal[im1][0]][0];
760 xm[ico][1]=fX[fIndLocal[im2][1]][1];
761 ym[ico][1]=fY[fIndLocal[im2][1]][1];
762
763 ixm[ico][0]=fIx[fIndLocal[im1][0]][0];
764 iym[ico][0]=fIy[fIndLocal[im1][0]][0];
765 ixm[ico][1]=fIx[fIndLocal[im2][1]][1];
766 iym[ico][1]=fIy[fIndLocal[im2][1]][1];
767 ico++;
768 }
769 }
770
771 Int_t nIco = ico;
772
773 fprintf(stderr,"nIco %d\n",nIco);
774 for (ico=0; ico<nIco; ico++) {
775 fprintf(stderr,"ico = %d\n",ico);
776 isec=fInput->Segmentation(0)->Sector(ixm[ico][0], iym[ico][0]);
777 dpx=fInput->Segmentation(0)->Dpx(isec)/2.;
778 dx=TMath::Abs(xm[ico][0]-xm[ico][1]);
779 isec=fInput->Segmentation(1)->Sector(ixm[ico][1], iym[ico][1]);
780 dpy=fInput->Segmentation(1)->Dpy(isec)/2.;
781 dy=TMath::Abs(ym[ico][0]-ym[ico][1]);
782
783 fprintf(stderr,"dx %f dpx %f dy %f dpy %f\n",dx,dpx,dy,dpy);
784 fprintf(stderr," X %f Y %f\n",xm[ico][1],ym[ico][0]);
785 if ((dx <= dpx) && (dy <= dpy)) {
786 fprintf(stderr,"ok\n");
787 Int_t cath;
788 AliMUONRawCluster cnew;
789 for (cath=0; cath<2; cath++) {
790 cnew.fX[cath]=Float_t(xm[ico][1]);
791 cnew.fY[cath]=Float_t(ym[ico][0]);
792 cnew.fMultiplicity[cath]=c->fMultiplicity[cath];
793 for (i=0; i<fMul[cath]; i++) {
794 cnew.fIndexMap[i][cath]=c->fIndexMap[i][cath];
795 fInput->Segmentation(cath)->SetPad(fIx[i][cath], fIy[i][cath]);
796 }
797 FillCluster(&cnew,cath);
798 }
799 cnew.fClusterType=cnew.PhysicsContribution();
800 AddRawCluster(cnew);
801 fNPeaks++;
802 }
803 }
804 delete [] xm;
805 delete [] ym;
806 delete [] ixm;
807 delete [] iym;
808 }
809}
810
811void AliMUONClusterFinderVS::FindLocalMaxima(AliMUONRawCluster* c)
812{
813// Find all local maxima of a cluster
814
815 AliMUONDigit* digt;
816
817 Int_t cath, cath1; // loops over cathodes
818 Int_t i; // loops over digits
819 Int_t j; // loops over cathodes
820//
821// Find local maxima
822//
823// counters for number of local maxima
824 fNLocal[0]=fNLocal[1]=0;
825// flags digits as local maximum
826 Bool_t isLocal[100][2];
827 for (i=0; i<100;i++) {
828 isLocal[i][0]=isLocal[i][1]=kFALSE;
829 }
830// number of next neighbours and arrays to store them
831 Int_t nn;
832 Int_t x[10], y[10];
833// loop over cathodes
834 for (cath=0; cath<2; cath++) {
835// loop over cluster digits
836 for (i=0; i<fMul[cath]; i++) {
837// get neighbours for that digit and assume that it is local maximum
838 fInput->Segmentation(cath)->Neighbours(fIx[i][cath], fIy[i][cath], &nn, x, y);
839 isLocal[i][cath]=kTRUE;
840 Int_t isec= fInput->Segmentation(cath)->Sector(fIx[i][cath], fIy[i][cath]);
841 Float_t a0 = fInput->Segmentation(cath)->Dpx(isec)*fInput->Segmentation(cath)->Dpy(isec);
842// loop over next neighbours, if at least one neighbour has higher charger assumption
843// digit is not local maximum
844 for (j=0; j<nn; j++) {
845 if (fHitMap[cath]->TestHit(x[j], y[j])==kEmpty) continue;
846 digt=(AliMUONDigit*) fHitMap[cath]->GetHit(x[j], y[j]);
847 isec=fInput->Segmentation(cath)->Sector(x[j], y[j]);
848 Float_t a1 = fInput->Segmentation(cath)->Dpx(isec)*fInput->Segmentation(cath)->Dpy(isec);
849 if (digt->fSignal/a1 > fQ[i][cath]/a0) {
850 isLocal[i][cath]=kFALSE;
851 break;
852//
853// handle special case of neighbouring pads with equal signal
854 } else if (digt->fSignal == fQ[i][cath]) {
855 if (fNLocal[cath]>0) {
856 for (Int_t k=0; k<fNLocal[cath]; k++) {
857 if (x[j]==fIx[fIndLocal[k][cath]][cath]
858 && y[j]==fIy[fIndLocal[k][cath]][cath])
859 {
860 isLocal[i][cath]=kFALSE;
861 }
862 } // loop over local maxima
863 } // are there already local maxima
864 } // same charge ?
865 } // loop over next neighbours
866 if (isLocal[i][cath]) {
867 fIndLocal[fNLocal[cath]][cath]=i;
868 fNLocal[cath]++;
869 }
870 } // loop over all digits
871 } // loop over cathodes
872
873 printf("\n Found %d %d %d %d local Maxima\n",
874 fNLocal[0], fNLocal[1], fMul[0], fMul[1]);
875 fprintf(stderr,"\n Cathode 1 local Maxima %d Multiplicite %d\n",fNLocal[0], fMul[0]);
876 fprintf(stderr," Cathode 2 local Maxima %d Multiplicite %d\n",fNLocal[1], fMul[1]);
877 Int_t ix, iy, isec;
878 Float_t dpx, dpy;
879
880
881 if (fNLocal[1]==2 && (fNLocal[0]==1 || fNLocal[0]==0)) {
882 Int_t iback=fNLocal[0];
883
884// Two local maxima on cathode 2 and one maximum on cathode 1
885// Look for local maxima considering up and down neighbours on the 1st cathode only
886//
887// Loop over cluster digits
888 cath=0;
889 cath1=1;
890
891 for (i=0; i<fMul[cath]; i++) {
892 isec=fInput->Segmentation(cath)->Sector(fIx[i][cath],fIy[i][cath]);
893 dpy=fInput->Segmentation(cath)->Dpy(isec);
894 dpx=fInput->Segmentation(cath)->Dpx(isec);
895 if (isLocal[i][cath]) continue;
896// Pad position should be consistent with position of local maxima on the opposite cathode
897 if ((TMath::Abs(fX[i][cath]-fX[fIndLocal[0][cath1]][cath1]) > dpx/2.) &&
898 (TMath::Abs(fX[i][cath]-fX[fIndLocal[1][cath1]][cath1]) > dpx/2.))
899 continue;
900
901// get neighbours for that digit and assume that it is local maximum
902 isLocal[i][cath]=kTRUE;
903// compare signal to that on the two neighbours on the left and on the right
904 fInput->Segmentation(cath)->GetPadI(fX[i][cath],fY[i][cath]+dpy,0,ix,iy);
905// iNN counts the number of neighbours with signal, it should be 1 or 2
906 Int_t iNN=0;
907 if (fHitMap[cath]->TestHit(ix, iy)!=kEmpty) {
908 iNN++;
909 digt=(AliMUONDigit*) fHitMap[cath]->GetHit(ix,iy);
910 if (digt->fSignal > fQ[i][cath]) isLocal[i][cath]=kFALSE;
911 }
912 fInput->Segmentation(cath)->GetPadI(fX[i][cath],fY[i][cath]-dpy,0,ix,iy);
913 if (fHitMap[cath]->TestHit(ix, iy)!=kEmpty) {
914 iNN++;
915 digt=(AliMUONDigit*) fHitMap[cath]->GetHit(ix,iy);
916 if (digt->fSignal > fQ[i][cath]) isLocal[i][cath]=kFALSE;
917 }
918 if (isLocal[i][cath] && iNN>0) {
919 fIndLocal[fNLocal[cath]][cath]=i;
920 fNLocal[cath]++;
921 }
922 } // loop over all digits
923// if one additional maximum has been found we are happy
924// if more maxima have been found restore the previous situation
925 fprintf(stderr,"\n New search gives %d local maxima for cathode 1 \n",fNLocal[0]);
926 fprintf(stderr," %d local maxima for cathode 2 \n",fNLocal[1]);
927 if (fNLocal[cath]>2) {
928 fNLocal[cath]=iback;
929 }
930
931 } // 1,2 local maxima
932
933 if (fNLocal[0]==2 && (fNLocal[1]==1 || fNLocal[1]==0)) {
934 Int_t iback=fNLocal[1];
935
936// Two local maxima on cathode 1 and one maximum on cathode 2
937// Look for local maxima considering left and right neighbours on the 2nd cathode only
938 cath=1;
939 Int_t cath1=0;
940
941
942//
943// Loop over cluster digits
944 for (i=0; i<fMul[cath]; i++) {
945 isec=fInput->Segmentation(cath)->Sector(fIx[i][cath],fIy[i][cath]);
946 dpx=fInput->Segmentation(cath)->Dpx(isec);
947 dpy=fInput->Segmentation(cath)->Dpy(isec);
948 if (isLocal[i][cath]) continue;
949// Pad position should be consistent with position of local maxima on the opposite cathode
950 if ((TMath::Abs(fY[i][cath]-fY[fIndLocal[0][cath1]][cath1]) > dpy/2.) &&
951 (TMath::Abs(fY[i][cath]-fY[fIndLocal[1][cath1]][cath1]) > dpy/2.))
952 continue;
953//
954// get neighbours for that digit and assume that it is local maximum
955 isLocal[i][cath]=kTRUE;
956// compare signal to that on the two neighbours on the left and on the right
957 fInput->Segmentation(cath)->GetPadI(fX[i][cath]+dpx,fY[i][cath],0,ix,iy);
958// iNN counts the number of neighbours with signal, it should be 1 or 2
959 Int_t iNN=0;
960 if (fHitMap[cath]->TestHit(ix, iy)!=kEmpty) {
961 iNN++;
962 digt=(AliMUONDigit*) fHitMap[cath]->GetHit(ix,iy);
963 if (digt->fSignal > fQ[i][cath]) isLocal[i][cath]=kFALSE;
964 }
965 fInput->Segmentation(cath)->GetPadI(fX[i][cath]-dpx,fY[i][cath],0,ix,iy);
966 if (fHitMap[cath]->TestHit(ix, iy)!=kEmpty) {
967 iNN++;
968 digt=(AliMUONDigit*) fHitMap[cath]->GetHit(ix,iy);
969 if (digt->fSignal > fQ[i][cath]) isLocal[i][cath]=kFALSE;
970 }
971 if (isLocal[i][cath] && iNN>0) {
972 fIndLocal[fNLocal[cath]][cath]=i;
973 fNLocal[cath]++;
974 }
975 } // loop over all digits
976// if one additional maximum has been found we are happy
977// if more maxima have been found restore the previous situation
978 fprintf(stderr,"\n New search gives %d local maxima for cathode 1 \n",fNLocal[0]);
979 fprintf(stderr,"\n %d local maxima for cathode 2 \n",fNLocal[1]);
980// printf("\n New search gives %d %d \n",fNLocal[0],fNLocal[1]);
981 if (fNLocal[cath]>2) {
982 fNLocal[cath]=iback;
983 }
984
985
986
987 } // 2,1 local maxima
988}
989
990
991void AliMUONClusterFinderVS::FillCluster(AliMUONRawCluster* c, Int_t flag, Int_t cath)
992{
993//
994// Completes cluster information starting from list of digits
995//
996 AliMUONDigit* dig;
997 Float_t x, y, z;
998 Int_t ix, iy;
999
1000 if (cath==1) {
1001 c->fPeakSignal[cath]=c->fPeakSignal[0];
1002 } else {
1003 c->fPeakSignal[cath]=0;
1004 }
1005
1006
1007 if (flag) {
1008 c->fX[cath]=0;
1009 c->fY[cath]=0;
1010 c->fQ[cath]=0;
1011 }
1012
1013// fprintf(stderr,"\n fPeakSignal %d\n",c->fPeakSignal[cath]);
1014 for (Int_t i=0; i<c->fMultiplicity[cath]; i++)
1015 {
1016 dig= fInput->Digit(cath,c->fIndexMap[i][cath]);
1017 ix=dig->fPadX+c->fOffsetMap[i][cath];
1018 iy=dig->fPadY;
1019 Int_t q=dig->fSignal;
1020 if (!flag) q=Int_t(q*c->fContMap[i][cath]);
1021// fprintf(stderr,"q %d c->fPeakSignal[ %d ] %d\n",q,cath,c->fPeakSignal[cath]);
1022 if (dig->fPhysics >= dig->fSignal) {
1023 c->fPhysicsMap[i]=2;
1024 } else if (dig->fPhysics == 0) {
1025 c->fPhysicsMap[i]=0;
1026 } else c->fPhysicsMap[i]=1;
1027//
1028//
1029// fprintf(stderr,"q %d c->fPeakSignal[cath] %d\n",q,c->fPeakSignal[cath]);
1030// peak signal and track list
1031 if (q>c->fPeakSignal[cath]) {
1032 c->fPeakSignal[cath]=q;
1033 c->fTracks[0]=dig->fHit;
1034 c->fTracks[1]=dig->fTracks[0];
1035 c->fTracks[2]=dig->fTracks[1];
1036// fprintf(stderr," c->fTracks[0] %d c->fTracks[1] %d\n",dig->fHit,dig->fTracks[0]);
1037 }
1038//
1039 if (flag) {
1040 fInput->Segmentation(cath)->GetPadC(ix, iy, x, y, z);
1041 c->fX[cath] += q*x;
1042 c->fY[cath] += q*y;
1043 c->fQ[cath] += q;
1044 }
1045 } // loop over digits
1046// fprintf(stderr," fin du cluster c\n");
1047
1048
1049 if (flag) {
1050 c->fX[cath]/=c->fQ[cath];
1051 c->fX[cath]=fInput->Segmentation(cath)->GetAnod(c->fX[cath]);
1052 c->fY[cath]/=c->fQ[cath];
1053//
1054// apply correction to the coordinate along the anode wire
1055//
1056 x=c->fX[cath];
1057 y=c->fY[cath];
1058 fInput->Segmentation(cath)->GetPadI(x, y, 0, ix, iy);
1059 fInput->Segmentation(cath)->GetPadC(ix, iy, x, y, z);
1060 Int_t isec=fInput->Segmentation(cath)->Sector(ix,iy);
1061 TF1* cogCorr = fInput->Segmentation(cath)->CorrFunc(isec-1);
1062
1063 if (cogCorr) {
1064 Float_t yOnPad=(c->fY[cath]-y)/fInput->Segmentation(cath)->Dpy(isec);
1065 c->fY[cath]=c->fY[cath]-cogCorr->Eval(yOnPad, 0, 0);
1066 }
1067 }
1068}
1069
1070void AliMUONClusterFinderVS::FillCluster(AliMUONRawCluster* c, Int_t cath)
1071{
1072//
1073// Completes cluster information starting from list of digits
1074//
1075 static Float_t dr0;
1076
1077 AliMUONDigit* dig;
1078
1079 if (cath==0) {
1080 dr0 = 10000;
1081 }
1082
1083 Float_t xpad, ypad, zpad;
1084 Float_t dx, dy, dr;
1085
1086 for (Int_t i=0; i<c->fMultiplicity[cath]; i++)
1087 {
1088 dig = fInput->Digit(cath,c->fIndexMap[i][cath]);
1089 fInput->Segmentation(cath)->
1090 GetPadC(dig->fPadX,dig->fPadY,xpad,ypad, zpad);
1091 fprintf(stderr,"x %f y %f cx %f cy %f\n",xpad,ypad,c->fX[0],c->fY[0]);
1092 dx = xpad - c->fX[0];
1093 dy = ypad - c->fY[0];
1094 dr = TMath::Sqrt(dx*dx+dy*dy);
1095
1096 if (dr < dr0) {
1097 dr0 = dr;
1098 fprintf(stderr," dr %f\n",dr);
1099 Int_t q=dig->fSignal;
1100 if (dig->fPhysics >= dig->fSignal) {
1101 c->fPhysicsMap[i]=2;
1102 } else if (dig->fPhysics == 0) {
1103 c->fPhysicsMap[i]=0;
1104 } else c->fPhysicsMap[i]=1;
1105 c->fPeakSignal[cath]=q;
1106 c->fTracks[0]=dig->fHit;
1107 c->fTracks[1]=dig->fTracks[0];
1108 c->fTracks[2]=dig->fTracks[1];
1109 fprintf(stderr," c->fTracks[0] %d c->fTracks[1] %d\n",dig->fHit,dig->fTracks[0]);
1110 }
1111//
1112 } // loop over digits
1113
1114// apply correction to the coordinate along the anode wire
1115 c->fX[cath]=fInput->Segmentation(cath)->GetAnod(c->fX[cath]);
1116}
1117
1118void AliMUONClusterFinderVS::FindCluster(Int_t i, Int_t j, Int_t cath, AliMUONRawCluster &c){
1119//
1120// Find clusterset
1121//
1122//
1123// Add i,j as element of the cluster
1124//
1125
1126 Int_t idx = fHitMap[cath]->GetHitIndex(i,j);
1127 AliMUONDigit* dig = (AliMUONDigit*) fHitMap[cath]->GetHit(i,j);
1128 Int_t q=dig->fSignal;
1129 Int_t theX=dig->fPadX;
1130 Int_t theY=dig->fPadY;
1131 if (q > TMath::Abs(c.fPeakSignal[0]) && q > TMath::Abs(c.fPeakSignal[1])) {
1132 c.fPeakSignal[cath]=q;
1133 c.fTracks[0]=dig->fHit;
1134 c.fTracks[1]=dig->fTracks[0];
1135 c.fTracks[2]=dig->fTracks[1];
1136 }
1137
1138//
1139// Make sure that list of digits is ordered
1140//
1141 Int_t mu=c.fMultiplicity[cath];
1142 c.fIndexMap[mu][cath]=idx;
1143
1144 if (dig->fPhysics >= dig->fSignal) {
1145 c.fPhysicsMap[mu]=2;
1146 } else if (dig->fPhysics == 0) {
1147 c.fPhysicsMap[mu]=0;
1148 } else c.fPhysicsMap[mu]=1;
1149 if (mu > 0) {
1150 for (Int_t ind=mu-1; ind>=0; ind--) {
1151 Int_t ist=(c.fIndexMap)[ind][cath];
1152 Int_t ql=fInput->Digit(cath, ist)->fSignal;
1153 Int_t ix=fInput->Digit(cath, ist)->fPadX;
1154 Int_t iy=fInput->Digit(cath, ist)->fPadY;
1155
1156 if (q>ql || (q==ql && theX > ix && theY < iy)) {
1157 c.fIndexMap[ind][cath]=idx;
1158 c.fIndexMap[ind+1][cath]=ist;
1159 } else {
1160 break;
1161 }
1162 }
1163 }
1164
1165 c.fMultiplicity[cath]++;
1166 if (c.fMultiplicity[cath] >= 50 ) {
1167 printf("FindCluster - multiplicity >50 %d \n",c.fMultiplicity[0]);
1168 c.fMultiplicity[cath]=49;
1169 }
1170
1171// Prepare center of gravity calculation
1172 Float_t x, y, z;
1173 fInput->Segmentation(cath)->GetPadC(i, j, x, y, z);
1174
1175 c.fX[cath] += q*x;
1176 c.fY[cath] += q*y;
1177 c.fQ[cath] += q;
1178// Flag hit as taken
1179 fHitMap[cath]->FlagHit(i,j);
1180//
1181// Now look recursively for all neighbours and pad hit on opposite cathode
1182//
1183// Loop over neighbours
1184 Int_t ix,iy;
1185 Int_t nn;
1186 Int_t xList[10], yList[10];
1187 fInput->Segmentation(cath)->Neighbours(i,j,&nn,xList,yList);
1188 for (Int_t in=0; in<nn; in++) {
1189 ix=xList[in];
1190 iy=yList[in];
1191 if (fHitMap[cath]->TestHit(ix,iy)==kUnused) FindCluster(ix, iy, cath, c);
1192 }
1193// Neighbours on opposite cathode
1194// Take into account that several pads can overlap with the present pad
1195 Float_t xmin, xmax, ymin, ymax, xc, yc;
1196 Int_t iop;
1197 Int_t isec=fInput->Segmentation(cath)->Sector(i,j);
1198 if (cath==0) {
1199 iop=1;
1200 xmin=x-fInput->Segmentation(cath)->Dpx(isec);
1201 xmax=x+fInput->Segmentation(cath)->Dpx(isec);
1202 xc=xmin+.001;
1203 while (xc < xmax) {
1204 xc+=fInput->Segmentation(iop)->Dpx(isec);
1205 fInput->Segmentation(iop)->GetPadI(xc,y,0,ix,iy);
1206 if (ix>=(fInput->Segmentation(iop)->Npx()) || (iy>=fInput->Segmentation(iop)->Npy())) continue;
1207 if (fHitMap[iop]->TestHit(ix,iy)==kUnused) FindCluster(ix, iy, iop, c);
1208 }
1209 } else {
1210 iop=0;
1211 ymin=y-fInput->Segmentation(cath)->Dpy(isec);
1212 ymax=y+fInput->Segmentation(cath)->Dpy(isec);
1213 yc=ymin+.001;
1214 while (yc < ymax) {
1215 yc+=fInput->Segmentation(iop)->Dpy(isec);
1216 fInput->Segmentation(iop)->GetPadI(x,yc,0,ix,iy);
1217 if (ix>=(fInput->Segmentation(iop)->Npx()) || (iy>=fInput->Segmentation(iop)->Npy())) continue;
1218 if (fHitMap[iop]->TestHit(ix,iy)==kUnused) FindCluster(ix, iy, iop, c);
1219 }
1220 }
1221}
1222
1223//_____________________________________________________________________________
1224
1225void AliMUONClusterFinderVS::FindRawClusters()
1226{
1227 //
1228 // MUON cluster finder from digits -- finds neighbours on both cathodes and
1229 // fills the tree with raw clusters
1230 //
1231
1232 if (!fInput->NDigits(0) && !fInput->NDigits(1)) return;
1233
1234 fHitMap[0] = new AliMUONHitMapA1(fInput->Segmentation(0), fInput->Digits(0));
1235 fHitMap[1] = new AliMUONHitMapA1(fInput->Segmentation(1), fInput->Digits(1));
1236
1237 AliMUONDigit *dig;
1238
1239 Int_t ndig, cath;
1240 Int_t nskip=0;
1241 Int_t ncls=0;
1242 fHitMap[0]->FillHits();
1243 fHitMap[1]->FillHits();
1244//
1245// Outer Loop over Cathodes
1246 for (cath=0; cath<2; cath++) {
1247 for (ndig=0; ndig<fInput->NDigits(cath); ndig++) {
1248 dig = fInput->Digit(cath, ndig);
1249 Int_t i=dig->fPadX;
1250 Int_t j=dig->fPadY;
1251 if (fHitMap[cath]->TestHit(i,j)==kUsed ||fHitMap[0]->TestHit(i,j)==kEmpty) {
1252 nskip++;
1253 continue;
1254 }
1255 fprintf(stderr,"\n CATHODE %d CLUSTER %d\n",cath,ncls);
1256 AliMUONRawCluster c;
1257 c.fMultiplicity[0]=0;
1258 c.fMultiplicity[1]=0;
1259 c.fPeakSignal[cath]=dig->fSignal;
1260 c.fTracks[0]=dig->fHit;
1261 c.fTracks[1]=dig->fTracks[0];
1262 c.fTracks[2]=dig->fTracks[1];
1263 // tag the beginning of cluster list in a raw cluster
1264 c.fNcluster[0]=-1;
1265
1266 FindCluster(i,j,cath,c);
1267
1268 // center of gravity
1269 c.fX[0] /= c.fQ[0];
1270 c.fX[0]=fInput->Segmentation(0)->GetAnod(c.fX[0]);
1271 c.fY[0] /= c.fQ[0];
1272 c.fX[1] /= c.fQ[1];
1273 c.fX[1]=fInput->Segmentation(0)->GetAnod(c.fX[1]);
1274 c.fY[1] /= c.fQ[1];
1275 fprintf(stderr,"\n Cathode 1 multiplicite %d X(CG) %f Y(CG) %f\n",c.fMultiplicity[0],c.fX[0],c.fY[0]);
1276 fprintf(stderr," Cathode 2 multiplicite %d X(CG) %f Y(CG) %f\n",c.fMultiplicity[1],c.fX[1],c.fY[1]);
1277
1278// Mathieson Fit
1279/*
1280 Bool_t fitted;
1281
1282 fitted=SingleMathiesonFit(&c, 0);
1283 c.fX[0]=fInput->Segmentation(0)->GetAnod(c.fX[0]);
1284 fitted=SingleMathiesonFit(&c, 1);
1285 c.fX[1]=fInput->Segmentation(1)->GetAnod(c.fX[1]);
1286*/
1287//
1288// Analyse cluster and decluster if necessary
1289//
1290 ncls++;
1291 c.fNcluster[1]=fNRawClusters;
1292 c.fClusterType=c.PhysicsContribution();
1293
1294 fNPeaks=0;
1295//
1296//
1297 Decluster(&c);
1298// AddRawCluster(c);
1299
1300//
1301// reset Cluster object
1302 { // begin local scope
1303 for (int k=0;k<c.fMultiplicity[0];k++) c.fIndexMap[k][0]=0;
1304 } // end local scope
1305
1306 { // begin local scope
1307 for (int k=0;k<c.fMultiplicity[1];k++) c.fIndexMap[k][1]=0;
1308 } // end local scope
1309
1310 c.fMultiplicity[0]=c.fMultiplicity[0]=0;
1311
1312
1313 } // end loop ndig
1314 } // end loop cathodes
1315 delete fHitMap[0];
1316 delete fHitMap[1];
1317}
1318
1319Float_t AliMUONClusterFinderVS::SingleMathiesonFit(AliMUONRawCluster *c, Int_t cath)
1320{
1321//
1322 AliMUONClusterInput& clusterInput = *(AliMUONClusterInput::Instance());
1323
1324 clusterInput.Fitter()->SetFCN(fcnS1);
1325 clusterInput.Fitter()->mninit(2,10,7);
1326 Double_t arglist[20];
1327 Int_t ierflag=0;
1328 arglist[0]=1;
1329// clusterInput.Fitter()->mnexcm("SET ERR",arglist,1,ierflag);
1330// Set starting values
1331 static Double_t vstart[2];
1332 vstart[0]=c->fX[1];
1333 vstart[1]=c->fY[0];
1334
1335
1336// lower and upper limits
1337 static Double_t lower[2], upper[2];
1338 Int_t ix,iy;
1339 fInput->Segmentation(cath)->GetPadI(c->fX[cath], c->fY[cath], 0, ix, iy);
1340 Int_t isec=fInput->Segmentation(cath)->Sector(ix, iy);
1341 lower[0]=vstart[0]-fInput->Segmentation(cath)->Dpx(isec)/2;
1342 lower[1]=vstart[1]-fInput->Segmentation(cath)->Dpy(isec)/2;
1343
1344 upper[0]=lower[0]+fInput->Segmentation(cath)->Dpx(isec);
1345 upper[1]=lower[1]+fInput->Segmentation(cath)->Dpy(isec);
1346
1347// step sizes
1348 static Double_t step[2]={0.0005, 0.0005};
1349
1350 clusterInput.Fitter()->mnparm(0,"x1",vstart[0],step[0],lower[0],upper[0],ierflag);
1351 clusterInput.Fitter()->mnparm(1,"y1",vstart[1],step[1],lower[1],upper[1],ierflag);
1352// ready for minimisation
1353 clusterInput.Fitter()->SetPrintLevel(1);
1354 clusterInput.Fitter()->mnexcm("SET OUT", arglist, 0, ierflag);
1355 arglist[0]= -1;
1356 arglist[1]= 0;
1357
1358 clusterInput.Fitter()->mnexcm("SET NOGR", arglist, 0, ierflag);
1359 clusterInput.Fitter()->mnexcm("MIGRAD", arglist, 0, ierflag);
1360 clusterInput.Fitter()->mnexcm("EXIT" , arglist, 0, ierflag);
1361 Double_t fmin, fedm, errdef;
1362 Int_t npari, nparx, istat;
1363
1364 clusterInput.Fitter()->mnstat(fmin, fedm, errdef, npari, nparx, istat);
1365 fFitStat=istat;
1366
1367// Print results
1368// Get fitted parameters
1369 Double_t xrec, yrec;
1370 TString chname;
1371 Double_t epxz, b1, b2;
1372 Int_t ierflg;
1373 clusterInput.Fitter()->mnpout(0, chname, xrec, epxz, b1, b2, ierflg);
1374 clusterInput.Fitter()->mnpout(1, chname, yrec, epxz, b1, b2, ierflg);
1375 fXFit[cath]=xrec;
1376 fYFit[cath]=yrec;
1377 return fmin;
1378}
1379
1380Float_t AliMUONClusterFinderVS::CombiSingleMathiesonFit(AliMUONRawCluster *c)
1381{
1382// Perform combined Mathieson fit on both cathode planes
1383//
1384 AliMUONClusterInput& clusterInput = *(AliMUONClusterInput::Instance());
1385 clusterInput.Fitter()->SetFCN(fcnCombiS1);
1386 clusterInput.Fitter()->mninit(2,10,7);
1387 Double_t arglist[20];
1388 Int_t ierflag=0;
1389 arglist[0]=1;
1390 static Double_t vstart[2];
1391 vstart[0]=fXInit[0];
1392 vstart[1]=fYInit[0];
1393
1394
1395// lower and upper limits
1396 static Double_t lower[2], upper[2];
1397 Int_t ix,iy,isec;
1398 fInput->Segmentation(0)->GetPadI(fXInit[0], fYInit[0], 0, ix, iy);
1399 isec=fInput->Segmentation(0)->Sector(ix, iy);
1400 Float_t dpy=fInput->Segmentation(0)->Dpy(isec)/2;
1401 fInput->Segmentation(1)->GetPadI(fXInit[0], fYInit[0], 0, ix, iy);
1402 isec=fInput->Segmentation(1)->Sector(ix, iy);
1403 Float_t dpx=fInput->Segmentation(1)->Dpx(isec)/2;
1404
1405
1406 lower[0]=vstart[0]-dpx;
1407 lower[1]=vstart[1]-dpy;
1408
1409 upper[0]=vstart[0]+dpx;
1410 upper[1]=vstart[1]+dpy;
1411
1412// step sizes
1413 static Double_t step[2]={0.00001, 0.0001};
1414
1415 clusterInput.Fitter()->mnparm(0,"x1",vstart[0],step[0],lower[0],upper[0],ierflag);
1416 clusterInput.Fitter()->mnparm(1,"y1",vstart[1],step[1],lower[1],upper[1],ierflag);
1417// ready for minimisation
1418 clusterInput.Fitter()->SetPrintLevel(1);
1419 clusterInput.Fitter()->mnexcm("SET OUT", arglist, 0, ierflag);
1420 arglist[0]= -1;
1421 arglist[1]= 0;
1422
1423 clusterInput.Fitter()->mnexcm("SET NOGR", arglist, 0, ierflag);
1424 clusterInput.Fitter()->mnexcm("MIGRAD", arglist, 0, ierflag);
1425 clusterInput.Fitter()->mnexcm("EXIT" , arglist, 0, ierflag);
1426 Double_t fmin, fedm, errdef;
1427 Int_t npari, nparx, istat;
1428
1429 clusterInput.Fitter()->mnstat(fmin, fedm, errdef, npari, nparx, istat);
1430 fFitStat=istat;
1431
1432// Print results
1433// Get fitted parameters
1434 Double_t xrec, yrec;
1435 TString chname;
1436 Double_t epxz, b1, b2;
1437 Int_t ierflg;
1438 clusterInput.Fitter()->mnpout(0, chname, xrec, epxz, b1, b2, ierflg);
1439 clusterInput.Fitter()->mnpout(1, chname, yrec, epxz, b1, b2, ierflg);
1440 fXFit[0]=xrec;
1441 fYFit[0]=yrec;
1442 return fmin;
1443}
1444
1445Bool_t AliMUONClusterFinderVS::DoubleMathiesonFit(AliMUONRawCluster *c, Int_t cath)
1446{
1447//
1448// Initialise global variables for fit
1449 AliMUONClusterInput& clusterInput = *(AliMUONClusterInput::Instance());
1450 clusterInput.Fitter()->SetFCN(fcnS2);
1451 clusterInput.Fitter()->mninit(5,10,7);
1452 Double_t arglist[20];
1453 Int_t ierflag=0;
1454 arglist[0]=1;
1455// Set starting values
1456 static Double_t vstart[5];
1457 vstart[0]=fX[fIndLocal[0][cath]][cath];
1458 vstart[1]=fY[fIndLocal[0][cath]][cath];
1459 vstart[2]=fX[fIndLocal[1][cath]][cath];
1460 vstart[3]=fY[fIndLocal[1][cath]][cath];
1461 vstart[4]=Float_t(fQ[fIndLocal[0][cath]][cath])/
1462 Float_t(fQ[fIndLocal[0][cath]][cath]+fQ[fIndLocal[1][cath]][cath]);
1463// lower and upper limits
1464 static Double_t lower[5], upper[5];
1465 Int_t isec=fInput->Segmentation(cath)->Sector(fIx[fIndLocal[0][cath]][cath], fIy[fIndLocal[0][cath]][cath]);
1466 lower[0]=vstart[0]-fInput->Segmentation(cath)->Dpx(isec);
1467 lower[1]=vstart[1]-fInput->Segmentation(cath)->Dpy(isec);
1468
1469 upper[0]=lower[0]+2.*fInput->Segmentation(cath)->Dpx(isec);
1470 upper[1]=lower[1]+2.*fInput->Segmentation(cath)->Dpy(isec);
1471
1472 isec=fInput->Segmentation(cath)->Sector(fIx[fIndLocal[1][cath]][cath], fIy[fIndLocal[1][cath]][cath]);
1473 lower[2]=vstart[2]-fInput->Segmentation(cath)->Dpx(isec)/2;
1474 lower[3]=vstart[3]-fInput->Segmentation(cath)->Dpy(isec)/2;
1475
1476 upper[2]=lower[2]+fInput->Segmentation(cath)->Dpx(isec);
1477 upper[3]=lower[3]+fInput->Segmentation(cath)->Dpy(isec);
1478
1479 lower[4]=0.;
1480 upper[4]=1.;
1481// step sizes
1482 static Double_t step[5]={0.0005, 0.0005, 0.0005, 0.0005, 0.0001};
1483
1484 clusterInput.Fitter()->mnparm(0,"x1",vstart[0],step[0],lower[0],upper[0],ierflag);
1485 clusterInput.Fitter()->mnparm(1,"y1",vstart[1],step[1],lower[1],upper[1],ierflag);
1486 clusterInput.Fitter()->mnparm(2,"x2",vstart[2],step[2],lower[2],upper[2],ierflag);
1487 clusterInput.Fitter()->mnparm(3,"y2",vstart[3],step[3],lower[3],upper[3],ierflag);
1488 clusterInput.Fitter()->mnparm(4,"a0",vstart[4],step[4],lower[4],upper[4],ierflag);
1489// ready for minimisation
1490 clusterInput.Fitter()->SetPrintLevel(-1);
1491 clusterInput.Fitter()->mnexcm("SET OUT", arglist, 0, ierflag);
1492 arglist[0]= -1;
1493 arglist[1]= 0;
1494
1495 clusterInput.Fitter()->mnexcm("SET NOGR", arglist, 0, ierflag);
1496 clusterInput.Fitter()->mnexcm("MIGRAD", arglist, 0, ierflag);
1497 clusterInput.Fitter()->mnexcm("EXIT" , arglist, 0, ierflag);
1498// Get fitted parameters
1499 Double_t xrec[2], yrec[2], qfrac;
1500 TString chname;
1501 Double_t epxz, b1, b2;
1502 Int_t ierflg;
1503 clusterInput.Fitter()->mnpout(0, chname, xrec[0], epxz, b1, b2, ierflg);
1504 clusterInput.Fitter()->mnpout(1, chname, yrec[0], epxz, b1, b2, ierflg);
1505 clusterInput.Fitter()->mnpout(2, chname, xrec[1], epxz, b1, b2, ierflg);
1506 clusterInput.Fitter()->mnpout(3, chname, yrec[1], epxz, b1, b2, ierflg);
1507 clusterInput.Fitter()->mnpout(4, chname, qfrac, epxz, b1, b2, ierflg);
1508
1509 Double_t fmin, fedm, errdef;
1510 Int_t npari, nparx, istat;
1511
1512 clusterInput.Fitter()->mnstat(fmin, fedm, errdef, npari, nparx, istat);
1513 fFitStat=istat;
1514 return kTRUE;
1515}
1516
1517Float_t AliMUONClusterFinderVS::CombiDoubleMathiesonFit(AliMUONRawCluster *c)
1518{
1519//
1520// Perform combined double Mathieson fit on both cathode planes
1521//
1522 AliMUONClusterInput& clusterInput = *(AliMUONClusterInput::Instance());
1523 clusterInput.Fitter()->SetFCN(fcnCombiS2);
1524 clusterInput.Fitter()->mninit(6,10,7);
1525 Double_t arglist[20];
1526 Int_t ierflag=0;
1527 arglist[0]=1;
1528// Set starting values
1529 static Double_t vstart[6];
1530 vstart[0]=fXInit[0];
1531 vstart[1]=fYInit[0];
1532 vstart[2]=fXInit[1];
1533 vstart[3]=fYInit[1];
1534 vstart[4]=fQrInit[0];
1535 vstart[5]=fQrInit[1];
1536// lower and upper limits
1537 static Double_t lower[6], upper[6];
1538 Int_t ix,iy,isec;
1539 Float_t dpx, dpy;
1540
1541 fInput->Segmentation(1)->GetPadI(fXInit[0], fYInit[0], 0, ix, iy);
1542 isec=fInput->Segmentation(1)->Sector(ix, iy);
1543 dpx=fInput->Segmentation(1)->Dpx(isec);
1544
1545 fInput->Segmentation(0)->GetPadI(fXInit[0], fYInit[0], 0, ix, iy);
1546 isec=fInput->Segmentation(0)->Sector(ix, iy);
1547 dpy=fInput->Segmentation(0)->Dpy(isec);
1548
1549 lower[0]=vstart[0]-dpx;
1550 lower[1]=vstart[1]-dpy;
1551 upper[0]=vstart[0]+dpx;
1552 upper[1]=vstart[1]+dpy;
1553
1554
1555 fInput->Segmentation(1)->GetPadI(fXInit[1], fYInit[1], 0, ix, iy);
1556 isec=fInput->Segmentation(1)->Sector(ix, iy);
1557 dpx=fInput->Segmentation(1)->Dpx(isec);
1558 fInput->Segmentation(0)->GetPadI(fXInit[1], fYInit[1], 0, ix, iy);
1559 isec=fInput->Segmentation(0)->Sector(ix, iy);
1560 dpy=fInput->Segmentation(0)->Dpy(isec);
1561
1562 lower[2]=vstart[2]-dpx;
1563 lower[3]=vstart[3]-dpy;
1564 upper[2]=vstart[2]+dpx;
1565 upper[3]=vstart[3]+dpy;
1566
1567
1568 lower[4]=0.;
1569 upper[4]=1.;
1570 lower[5]=0.;
1571 upper[5]=1.;
1572
1573// step sizes
1574 static Double_t step[6]={0.0005, 0.0005, 0.0005, 0.0005, 0.001, 0.001};
1575
1576 clusterInput.Fitter()->mnparm(0,"x1",vstart[0],step[0],lower[0],upper[0],ierflag);
1577 clusterInput.Fitter()->mnparm(1,"y1",vstart[1],step[1],lower[1],upper[1],ierflag);
1578 clusterInput.Fitter()->mnparm(2,"x2",vstart[2],step[2],lower[2],upper[2],ierflag);
1579 clusterInput.Fitter()->mnparm(3,"y2",vstart[3],step[3],lower[3],upper[3],ierflag);
1580 clusterInput.Fitter()->mnparm(4,"a0",vstart[4],step[4],lower[4],upper[4],ierflag);
1581 clusterInput.Fitter()->mnparm(5,"a1",vstart[5],step[5],lower[5],upper[5],ierflag);
1582// ready for minimisation
1583 clusterInput.Fitter()->SetPrintLevel(-1);
1584 clusterInput.Fitter()->mnexcm("SET OUT", arglist, 0, ierflag);
1585 arglist[0]= -1;
1586 arglist[1]= 0;
1587
1588 clusterInput.Fitter()->mnexcm("SET NOGR", arglist, 0, ierflag);
1589 clusterInput.Fitter()->mnexcm("MIGRAD", arglist, 0, ierflag);
1590 clusterInput.Fitter()->mnexcm("EXIT" , arglist, 0, ierflag);
1591// Get fitted parameters
1592 TString chname;
1593 Double_t epxz, b1, b2;
1594 Int_t ierflg;
1595 clusterInput.Fitter()->mnpout(0, chname, fXFit[0], epxz, b1, b2, ierflg);
1596 clusterInput.Fitter()->mnpout(1, chname, fYFit[0], epxz, b1, b2, ierflg);
1597 clusterInput.Fitter()->mnpout(2, chname, fXFit[1], epxz, b1, b2, ierflg);
1598 clusterInput.Fitter()->mnpout(3, chname, fYFit[1], epxz, b1, b2, ierflg);
1599 clusterInput.Fitter()->mnpout(4, chname, fQrFit[0], epxz, b1, b2, ierflg);
1600 clusterInput.Fitter()->mnpout(5, chname, fQrFit[1], epxz, b1, b2, ierflg);
1601
1602 Double_t fmin, fedm, errdef;
1603 Int_t npari, nparx, istat;
1604
1605 clusterInput.Fitter()->mnstat(fmin, fedm, errdef, npari, nparx, istat);
1606 fFitStat=istat;
1607
1608 fChi2[0]=fmin;
1609 fChi2[1]=fmin;
1610 return fmin;
1611}
1612
1613void AliMUONClusterFinderVS::Split(AliMUONRawCluster* c)
1614{
1615//
1616// One cluster for each maximum
1617//
1618 Int_t i, j, cath;
1619 AliMUONClusterInput& clusterInput = *(AliMUONClusterInput::Instance());
1620 for (j=0; j<2; j++) {
1621 AliMUONRawCluster cnew;
1622 for (cath=0; cath<2; cath++) {
1623 cnew.fChi2[cath]=fChi2[0];
1624
1625 if (fNPeaks == 0) {
1626 cnew.fNcluster[0]=-1;
1627 cnew.fNcluster[1]=fNRawClusters;
1628 } else {
1629 cnew.fNcluster[0]=fNPeaks;
1630 cnew.fNcluster[1]=0;
1631 }
1632 cnew.fMultiplicity[cath]=0;
1633 cnew.fX[cath]=Float_t(fXFit[j]);
1634 cnew.fY[cath]=Float_t(fYFit[j]);
1635 if (j==0) {
1636 cnew.fQ[cath]=Int_t(clusterInput.TotalCharge(cath)*fQrFit[cath]);
1637 } else {
1638 cnew.fQ[cath]=Int_t(clusterInput.TotalCharge(cath)*(1-fQrFit[cath]));
1639 }
1640 fInput->Segmentation(cath)->SetHit(fXFit[j],fYFit[j],0);
1641 for (i=0; i<fMul[cath]; i++) {
1642 cnew.fIndexMap[cnew.fMultiplicity[cath]][cath]=
1643 c->fIndexMap[i][cath];
1644 fInput->Segmentation(cath)->SetPad(fIx[i][cath], fIy[i][cath]);
1645 Float_t q1=fInput->Response()->IntXY(fInput->Segmentation(cath));
1646 cnew.fContMap[i][cath]
1647 =(q1*Float_t(cnew.fQ[cath]))/Float_t(fQ[i][cath]);
1648 cnew.fMultiplicity[cath]++;
1649// printf(" fXFIT %f fYFIT %f Multiplicite %d\n",cnew.fX[cath],cnew.fY[cath],cnew.fMultiplicity[cath]);
1650 }
1651 FillCluster(&cnew,0,cath);
1652 } // cathode loop
1653
1654 cnew.fClusterType=cnew.PhysicsContribution();
1655 if (cnew.fQ[0]>0 && cnew.fQ[1]>0) AddRawCluster(cnew);
1656 fNPeaks++;
1657 }
1658}
1659
1660
1661//
1662// Minimisation functions
1663// Single Mathieson
1664void fcnS1(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag)
1665{
1666 AliMUONClusterInput& clusterInput = *(AliMUONClusterInput::Instance());
1667 Int_t i;
1668 Float_t delta;
1669 Float_t chisq=0;
1670 Float_t qcont=0;
1671 Float_t qtot=0;
1672
1673 for (i=0; i<clusterInput.Nmul(0); i++) {
1674 Float_t q0=clusterInput.Charge(i,0);
1675 Float_t q1=clusterInput.DiscrChargeS1(i,par);
1676 delta=(q0-q1)/q0;
1677 chisq+=delta*delta;
1678 qcont+=q1;
1679 qtot+=q0;
1680 }
1681 f=chisq;
1682}
1683
1684void fcnCombiS1(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag)
1685{
1686 AliMUONClusterInput& clusterInput = *(AliMUONClusterInput::Instance());
1687 Int_t i, cath;
1688 Float_t delta;
1689 Float_t chisq=0;
1690 Float_t qcont=0;
1691 Float_t qtot=0;
1692 // Float_t chi2temp=0;
1693
1694 for (cath=0; cath<2; cath++) {
1695// chisq=0;
1696 for (i=0; i<clusterInput.Nmul(cath); i++) {
1697 Float_t q0=clusterInput.Charge(i,cath);
1698 Float_t q1=clusterInput.DiscrChargeCombiS1(i,par,cath);
1699 // delta=(q0-q1);
1700 delta=(q0-q1)/q0;
1701 chisq+=delta*delta;
1702 qcont+=q1;
1703 qtot+=q0;
1704 }
1705// if (cath == 0) chi2temp=chisq/clusterInput.Nbins[cath];
1706 }
1707// chisq = chisq/clusterInput.Nbins[1]+chi2temp;
1708 f=chisq;
1709}
1710
1711// Double Mathieson
1712void fcnS2(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag)
1713{
1714 AliMUONClusterInput& clusterInput = *(AliMUONClusterInput::Instance());
1715 Int_t i;
1716 Float_t delta;
1717 Float_t chisq=0;
1718 Float_t qcont=0;
1719 Float_t qtot=0;
1720
1721 for (i=0; i<clusterInput.Nmul(0); i++) {
1722
1723 Float_t q0=clusterInput.Charge(i,0);
1724 Float_t q1=clusterInput.DiscrChargeS2(i,par);
1725 delta=(q0-q1)/q0;
1726 chisq+=delta*delta;
1727 qcont+=q1;
1728 qtot+=q0;
1729 }
1730// chisq=chisq+=(qtot-qcont)*(qtot-qcont)*0.5;
1731 f=chisq;
1732}
1733
1734// Double Mathieson
1735void fcnCombiS2(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag)
1736{
1737 AliMUONClusterInput& clusterInput = *(AliMUONClusterInput::Instance());
1738 Int_t i, cath;
1739 Float_t delta;
1740 Float_t chisq=0;
1741 Float_t qcont=0;
1742 Float_t qtot=0;
1743 // Float_t chi2temp=0;
1744
1745 for (cath=0; cath<2; cath++) {
1746// chisq=0;
1747 for (i=0; i<clusterInput.Nmul(cath); i++) {
1748 Float_t q0=clusterInput.Charge(i,cath);
1749 Float_t q1=clusterInput.DiscrChargeCombiS2(i,par,cath);
1750 // delta=(q0-q1);
1751 delta=(q0-q1)/q0;
1752 chisq+=delta*delta;
1753 qcont+=q1;
1754 qtot+=q0;
1755 }
1756// if (cath == 0) chi2temp=chisq/clusterInput.Nbins[cath];
1757 }
1758// chisq = chisq/clusterInput.Nbins[1]+chi2temp;
1759 f=chisq;
1760}
1761
1762void AliMUONClusterFinderVS::AddRawCluster(const AliMUONRawCluster c)
1763{
1764 //
1765 // Add a raw cluster copy to the list
1766 //
1767 AliMUON *pMUON=(AliMUON*)gAlice->GetModule("MUON");
1768 pMUON->AddRawCluster(fInput->Chamber(),c);
1769 fNRawClusters++;
1770 fprintf(stderr,"\nfNRawClusters %d\n",fNRawClusters);
1771}
1772
1773Bool_t AliMUONClusterFinderVS::TestTrack(Int_t t) {
1774 if (fTrack[0]==-1 || fTrack[1]==-1) {
1775 return kTRUE;
1776 } else if (t==fTrack[0] || t==fTrack[1]) {
1777 return kTRUE;
1778 } else {
1779 return kFALSE;
1780 }
1781}
1782
1783AliMUONClusterFinderVS& AliMUONClusterFinderVS
1784::operator = (const AliMUONClusterFinderVS& rhs)
1785{
1786// Dummy assignment operator
1787 return *this;
1788}
1789
1790
1791
1792
1793
1794
1795
1796
1797