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