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