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