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