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