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