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