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