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