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