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