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