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