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