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