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