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