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