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