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