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