]> git.uio.no Git - u/mrichter/AliRoot.git/blame - MUON/AliMUONClusterFinder.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Introduction of the reference to Copyright and cvs Id
[u/mrichter/AliRoot.git] / MUON / AliMUONClusterFinder.cxx
CommitLineData
4c039060 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/*
17$Log$
18*/
19
a897a37a 20#include "AliMUONClusterFinder.h"
21#include "TTree.h"
22#include "AliRun.h"
23#include <TCanvas.h>
24#include <TH1.h>
25#include <TPad.h>
26#include <TGraph.h>
27#include <TPostScript.h>
28#include <TMinuit.h>
29
e3a4d40e 30//----------------------------------------------------------
a897a37a 31static AliMUONsegmentation* gSegmentation;
32static AliMUONresponse* gResponse;
33static Int_t gix[500];
34static Int_t giy[500];
35static Float_t gCharge[500];
36static Int_t gNbins;
37static Int_t gFirst=kTRUE;
38static TMinuit *gMyMinuit ;
39void fcn(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag);
40static Int_t gChargeTot;
41
a897a37a 42//----------------------------------------------------------
e3a4d40e 43
a897a37a 44ClassImp(AliMUONClusterFinder)
45
46 AliMUONClusterFinder::AliMUONClusterFinder
47(AliMUONsegmentation *segmentation, AliMUONresponse *response,
48 TClonesArray *digits, Int_t chamber)
49{
50 fSegmentation=segmentation;
51 fResponse=response;
52
53 fDigits=digits;
54 fNdigits = fDigits->GetEntriesFast();
55 fChamber=chamber;
56 fRawClusters=new TClonesArray("AliMUONRawCluster",10000);
57 fNRawClusters=0;
58 fCogCorr = 0;
59 SetNperMax();
60 SetClusterSize();
61 SetDeclusterFlag();
62 fNPeaks=-1;
63}
64
65 AliMUONClusterFinder::AliMUONClusterFinder()
66{
67 fSegmentation=0;
68 fResponse=0;
69
70 fDigits=0;
71 fNdigits = 0;
72 fChamber=-1;
73 fRawClusters=new TClonesArray("AliMUONRawCluster",10000);
74 fNRawClusters=0;
75 fHitMap = 0;
76 fCogCorr = 0;
77 SetNperMax();
78 SetClusterSize();
79 SetDeclusterFlag();
80 fNPeaks=-1;
81}
82
83void AliMUONClusterFinder::AddRawCluster(const AliMUONRawCluster c)
84{
85 //
86 // Add a raw cluster copy to the list
87 //
88 AliMUON *MUON=(AliMUON*)gAlice->GetModule("MUON");
89 MUON->AddRawCluster(fChamber,c);
90 fNRawClusters++;
91}
92
93
94
95void AliMUONClusterFinder::Decluster(AliMUONRawCluster *cluster)
96{
97// AliMUONdigit *dig;
98// Int_t q;
99
100
101 Int_t mul = cluster->fMultiplicity;
102// printf("Decluster - multiplicity %d \n",mul);
103
104 if (mul == 1 || mul ==2) {
105// printf("\n Nothing special for 1- and 2-clusters \n");
106//
107// Nothing special for 1- and 2-clusters
108 if (fNPeaks != 0) {
109 cluster->fNcluster[0]=fNPeaks;
110 cluster->fNcluster[1]=0;
111 }
112 AddRawCluster(*cluster);
113 fNPeaks++;
114 } else if (mul ==3) {
115//
116// 3-cluster, check topology
117// printf("\n 3-cluster, check topology \n");
118 if (fDeclusterFlag) {
119 if (Centered(cluster)) {
120 // ok, cluster is centered
121 // printf("\n ok, cluster is centered \n");
122 } else {
123 // cluster is not centered, split into 2+1
124 // printf("\n cluster is not centered, split into 2+1 \n");
125 }
126 } else {
127 if (fNPeaks != 0) {
128 cluster->fNcluster[0]=fNPeaks;
129 cluster->fNcluster[1]=0;
130 }
131 AddRawCluster(*cluster);
132 fNPeaks++;
133 }
134 } else {
135 //if (mul < 12) {
136 // printf("Decluster - multiplicity > 45 %d \n",mul);
137 //printf("Decluster - multiplicity < 25 %d \n",mul);
138//
139// 4-and more-pad clusters
140//
141 if (mul <= fClusterSize) {
142 if (fDeclusterFlag) {
143 SplitByLocalMaxima(cluster);
144 } else {
145 if (fNPeaks != 0) {
146 cluster->fNcluster[0]=fNPeaks;
147 cluster->fNcluster[1]=0;
148 }
149 AddRawCluster(*cluster);
150 fNPeaks++;
151 }
152 }
153 //}
154 } // multiplicity
155}
156
157
158Bool_t AliMUONClusterFinder::Centered(AliMUONRawCluster *cluster)
159{
160 AliMUONdigit* dig;
161 dig= (AliMUONdigit*)fDigits->UncheckedAt(cluster->fIndexMap[0]);
162 Int_t ix=dig->fPadX;
163 Int_t iy=dig->fPadY;
164 Int_t nn;
165 Int_t X[kMaxNeighbours], Y[kMaxNeighbours], XN[kMaxNeighbours], YN[kMaxNeighbours];
166
167 fSegmentation->Neighbours(ix,iy,&nn,X,Y);
168 Int_t nd=0;
169 for (Int_t i=0; i<nn; i++) {
170 if (fHitMap->TestHit(X[i],Y[i]) == used) {
171 XN[nd]=X[i];
172 YN[nd]=Y[i];
173 nd++;
174 }
175 }
176 if (nd==2) {
177//
178// cluster is centered !
179 if (fNPeaks != 0) {
180 cluster->fNcluster[0]=fNPeaks;
181 cluster->fNcluster[1]=0;
182 }
183 AddRawCluster(*cluster);
184 fNPeaks++;
185 return kTRUE;
186 } else if (nd ==1) {
187//
188// Highest signal on an edge, split cluster into 2+1
189//
190// who is the neighbour ?
191 Int_t nind=fHitMap->GetHitIndex(XN[0], YN[0]);
192 Int_t i1= (nind==cluster->fIndexMap[1]) ? 1:2;
193 Int_t i2= (nind==cluster->fIndexMap[1]) ? 2:1;
194//
195// 2-cluster
196 AliMUONRawCluster cnew;
197 if (fNPeaks == 0) {
198 cnew.fNcluster[0]=-1;
199 cnew.fNcluster[1]=fNRawClusters;
200 } else {
201 cnew.fNcluster[0]=fNPeaks;
202 cnew.fNcluster[1]=0;
203 }
204 cnew.fMultiplicity=2;
205 cnew.fIndexMap[0]=cluster->fIndexMap[0];
206 cnew.fIndexMap[1]=cluster->fIndexMap[i1];
207 FillCluster(&cnew);
208 cnew.fClusterType=cnew.PhysicsContribution();
209 AddRawCluster(cnew);
210 fNPeaks++;
211//
212// 1-cluster
213 cluster->fMultiplicity=1;
214 cluster->fIndexMap[0]=cluster->fIndexMap[i2];
215 cluster->fIndexMap[1]=0;
216 cluster->fIndexMap[2]=0;
217 FillCluster(cluster);
218 if (fNPeaks != 0) {
219 cluster->fNcluster[0]=fNPeaks;
220 cluster->fNcluster[1]=0;
221 }
222 cluster->fClusterType=cluster->PhysicsContribution();
223 AddRawCluster(*cluster);
224 fNPeaks++;
225 return kFALSE;
226 } else {
227 printf("\n Completely screwed up %d !! \n",nd);
228
229 }
230
231 return kFALSE;
232}
233void AliMUONClusterFinder::SplitByLocalMaxima(AliMUONRawCluster *c)
234{
235 AliMUONdigit* dig[100], *digt;
236 Int_t ix[100], iy[100], q[100];
237 Float_t x[100], y[100];
238 Int_t i; // loops over digits
239 Int_t j; // loops over local maxima
240 // Float_t xPeak[2];
241 // Float_t yPeak[2];
242 // Int_t threshold=500;
243 Int_t mul=c->fMultiplicity;
244//
245// dump digit information into arrays
246//
247 for (i=0; i<mul; i++)
248 {
249 dig[i]= (AliMUONdigit*)fDigits->UncheckedAt(c->fIndexMap[i]);
250 ix[i]= dig[i]->fPadX;
251 iy[i]= dig[i]->fPadY;
252 q[i] = dig[i]->fSignal;
253 fSegmentation->GetPadCxy(ix[i], iy[i], x[i], y[i]);
254 }
255//
256// Find local maxima
257//
258 Bool_t IsLocal[100];
259 Int_t NLocal=0;
260 Int_t AssocPeak[100];
261 Int_t IndLocal[100];
262 Int_t nn;
263 Int_t X[kMaxNeighbours], Y[kMaxNeighbours];
264 for (i=0; i<mul; i++) {
265 fSegmentation->Neighbours(ix[i], iy[i], &nn, X, Y);
266 IsLocal[i]=kTRUE;
267 for (j=0; j<nn; j++) {
e3a4d40e 268 if (fHitMap->TestHit(X[j],Y[j])==empty) continue;
a897a37a 269 digt=(AliMUONdigit*) fHitMap->GetHit(X[j], Y[j]);
270 if (digt->fSignal > q[i]) {
271 IsLocal[i]=kFALSE;
272 break;
273//
274// handle special case of neighbouring pads with equal signal
275 } else if (digt->fSignal == q[i]) {
276 if (NLocal >0) {
277 for (Int_t k=0; k<NLocal; k++) {
278 if (X[j]==ix[IndLocal[k]] && Y[j]==iy[IndLocal[k]]){
279 IsLocal[i]=kFALSE;
280 }
281 }
282 }
283 }
284 } // loop over next neighbours
285 // Maxima should not be on the edge
286 if (IsLocal[i]) {
287 IndLocal[NLocal]=i;
288 NLocal++;
289 }
290 } // loop over all digits
291// printf("Found %d local Maxima",NLocal);
292//
293// If only one local maximum found but multiplicity is high
294// take global maximum from the list of digits.
295 if (NLocal==1 && mul>12) {
296 Int_t nnew=0;
297 for (i=0; i<mul; i++) {
298 if (!IsLocal[i]) {
299 IndLocal[NLocal]=i;
300 IsLocal[i]=kTRUE;
301 NLocal++;
302 nnew++;
303 }
304 if (nnew==1) break;
305 }
306 }
307
308// If number of local maxima is 2 try to fit a double gaussian
309 if (NLocal==2) {
310//
311// Initialise global variables for fit
312 gFirst=1;
313 gSegmentation=fSegmentation;
314 gResponse =fResponse;
315 gNbins=mul;
316
317 for (i=0; i<mul; i++) {
318 gix[i]=ix[i];
319 giy[i]=iy[i];
320 gCharge[i]=Float_t(q[i]);
321 }
322//
323 if (gFirst) {
324 gFirst=kFALSE;
325 gMyMinuit = new TMinuit(5);
326 }
327 gMyMinuit->SetFCN(fcn);
328 gMyMinuit->mninit(5,10,7);
329 Double_t arglist[20];
330 Int_t ierflag=0;
331 arglist[0]=1;
332// gMyMinuit->mnexcm("SET ERR",arglist,1,ierflag);
333// Set starting values
334 static Double_t vstart[5];
335 vstart[0]=x[IndLocal[0]];
336 vstart[1]=y[IndLocal[0]];
337 vstart[2]=x[IndLocal[1]];
338 vstart[3]=y[IndLocal[1]];
339 vstart[4]=Float_t(q[IndLocal[0]])/
340 Float_t(q[IndLocal[0]]+q[IndLocal[1]]);
341// lower and upper limits
342 static Double_t lower[5], upper[5];
343 Int_t isec=fSegmentation->Sector(ix[IndLocal[0]], iy[IndLocal[0]]);
344 lower[0]=vstart[0]-fSegmentation->Dpx(isec)/2;
345 lower[1]=vstart[1]-fSegmentation->Dpy(isec)/2;
346// lower[1]=vstart[1];
347
348 upper[0]=lower[0]+fSegmentation->Dpx(isec);
349 upper[1]=lower[1]+fSegmentation->Dpy(isec);
350// upper[1]=vstart[1];
351
352 isec=fSegmentation->Sector(ix[IndLocal[1]], iy[IndLocal[1]]);
353 lower[2]=vstart[2]-fSegmentation->Dpx(isec)/2;
354 lower[3]=vstart[3]-fSegmentation->Dpy(isec)/2;
355// lower[3]=vstart[3];
356
357 upper[2]=lower[2]+fSegmentation->Dpx(isec);
358 upper[3]=lower[3]+fSegmentation->Dpy(isec);
359// upper[3]=vstart[3];
360
361 lower[4]=0.;
362 upper[4]=1.;
363// step sizes
364 static Double_t step[5]={0.005, 0.03, 0.005, 0.03, 0.01};
365
366 gMyMinuit->mnparm(0,"x1",vstart[0],step[0],lower[0],upper[0],ierflag);
367 gMyMinuit->mnparm(1,"y1",vstart[1],step[1],lower[1],upper[1],ierflag);
368 gMyMinuit->mnparm(2,"x2",vstart[2],step[2],lower[2],upper[2],ierflag);
369 gMyMinuit->mnparm(3,"y2",vstart[3],step[3],lower[3],upper[3],ierflag);
370 gMyMinuit->mnparm(4,"a0",vstart[4],step[4],lower[4],upper[4],ierflag);
371// ready for minimisation
372 gMyMinuit->SetPrintLevel(-1);
373 gMyMinuit->mnexcm("SET OUT", arglist, 0, ierflag);
374 arglist[0]= -1;
375 arglist[1]= 0;
376
377 gMyMinuit->mnexcm("SET NOGR", arglist, 0, ierflag);
378 gMyMinuit->mnexcm("SCAN", arglist, 0, ierflag);
379 gMyMinuit->mnexcm("EXIT" , arglist, 0, ierflag);
380// Print results
381// Double_t amin,edm,errdef;
382// Int_t nvpar,nparx,icstat;
383// gMyMinuit->mnstat(amin,edm,errdef,nvpar,nparx,icstat);
384// gMyMinuit->mnprin(3,amin);
385// Get fitted parameters
386
387 Double_t xrec[2], yrec[2], qfrac;
388 TString chname;
389 Double_t epxz, b1, b2;
390 Int_t ierflg;
391 gMyMinuit->mnpout(0, chname, xrec[0], epxz, b1, b2, ierflg);
392 gMyMinuit->mnpout(1, chname, yrec[0], epxz, b1, b2, ierflg);
393 gMyMinuit->mnpout(2, chname, xrec[1], epxz, b1, b2, ierflg);
394 gMyMinuit->mnpout(3, chname, yrec[1], epxz, b1, b2, ierflg);
395 gMyMinuit->mnpout(4, chname, qfrac, epxz, b1, b2, ierflg);
396 printf("\n %f %f %f %f %f\n", xrec[0], yrec[0], xrec[1], yrec[1],qfrac);
397// delete gMyMinuit;
398
399
400 //
401 // One cluster for each maximum
402 //
403 for (j=0; j<2; j++) {
404 AliMUONRawCluster cnew;
405 if (fNPeaks == 0) {
406 cnew.fNcluster[0]=-1;
407 cnew.fNcluster[1]=fNRawClusters;
408 } else {
409 cnew.fNcluster[0]=fNPeaks;
410 cnew.fNcluster[1]=0;
411 }
412 cnew.fMultiplicity=0;
413 cnew.fX=Float_t(xrec[j]);
414 cnew.fY=Float_t(yrec[j]);
415 if (j==0) {
a6f39961 416 cnew.fQ=Int_t(gChargeTot*qfrac);
a897a37a 417 } else {
a6f39961 418 cnew.fQ=Int_t(gChargeTot*(1-qfrac));
a897a37a 419 }
420 gSegmentation->SetHit(xrec[j],yrec[j]);
421 for (i=0; i<mul; i++) {
422 cnew.fIndexMap[cnew.fMultiplicity]=c->fIndexMap[i];
423 gSegmentation->SetPad(gix[i], giy[i]);
424 Float_t q1=gResponse->IntXY(gSegmentation);
425 cnew.fContMap[cnew.fMultiplicity]=Float_t(q[i])/(q1*cnew.fQ);
426 cnew.fMultiplicity++;
427 }
428 FillCluster(&cnew,0);
429 //printf("\n x,y %f %f ", cnew.fX, cnew.fY);
430 cnew.fClusterType=cnew.PhysicsContribution();
431 AddRawCluster(cnew);
432 fNPeaks++;
433 }
434 }
435
a6f39961 436 Bool_t fitted=kTRUE;
a897a37a 437
438 if (NLocal !=2 || !fitted) {
439 // Check if enough local clusters have been found,
440 // if not add global maxima to the list
441 //
442 Int_t nPerMax;
443 if (NLocal!=0) {
444 nPerMax=mul/NLocal;
445 } else {
446 printf("\n Warning, no local maximum found \n");
447 nPerMax=fNperMax+1;
448 }
449
450 if (nPerMax > fNperMax) {
451 Int_t nGlob=mul/fNperMax-NLocal+1;
452 if (nGlob > 0) {
453 Int_t nnew=0;
454 for (i=0; i<mul; i++) {
455 if (!IsLocal[i]) {
456 IndLocal[NLocal]=i;
457 IsLocal[i]=kTRUE;
458 NLocal++;
459 nnew++;
460 }
461 if (nnew==nGlob) break;
462 }
463 }
464 }
465 //
466 // Associate hits to peaks
467 //
468 for (i=0; i<mul; i++) {
469 Float_t dmin=1.E10;
470 Float_t qmax=0;
471 if (IsLocal[i]) continue;
472 for (j=0; j<NLocal; j++) {
473 Int_t il=IndLocal[j];
474 Float_t d=TMath::Sqrt((x[i]-x[il])*(x[i]-x[il])
475 +(y[i]-y[il])*(y[i]-y[il]));
476 Float_t ql=q[il];
477 //
478 // Select nearest peak
479 //
480 if (d<dmin) {
481 dmin=d;
482 qmax=ql;
483 AssocPeak[i]=j;
484 } else if (d==dmin) {
485 //
486 // If more than one take highest peak
487 //
488 if (ql>qmax) {
489 dmin=d;
490 qmax=ql;
491 AssocPeak[i]=j;
492 }
493 }
494 }
495 }
496
497
498 //
499 // One cluster for each maximum
500 //
501 for (j=0; j<NLocal; j++) {
502 AliMUONRawCluster cnew;
503 if (fNPeaks == 0) {
504 cnew.fNcluster[0]=-1;
505 cnew.fNcluster[1]=fNRawClusters;
506 } else {
507 cnew.fNcluster[0]=fNPeaks;
508 cnew.fNcluster[1]=0;
509 }
510 cnew.fIndexMap[0]=c->fIndexMap[IndLocal[j]];
511 cnew.fMultiplicity=1;
512 for (i=0; i<mul; i++) {
513 if (IsLocal[i]) continue;
514 if (AssocPeak[i]==j) {
515 cnew.fIndexMap[cnew.fMultiplicity]=c->fIndexMap[i];
516 cnew.fMultiplicity++;
517 }
518 }
519 FillCluster(&cnew);
520 cnew.fClusterType=cnew.PhysicsContribution();
521 AddRawCluster(cnew);
522 fNPeaks++;
523 }
524 }
525
526}
527
528
529void AliMUONClusterFinder::FillCluster(AliMUONRawCluster* c, Int_t flag)
530{
531//
532// Completes cluster information starting from list of digits
533//
534 AliMUONdigit* dig;
535 Float_t x, y;
536 Int_t ix, iy;
537 Float_t frac=0;
538
539 c->fPeakSignal=0;
540 if (flag) {
541 c->fX=0;
542 c->fY=0;
543 c->fQ=0;
544 }
545 //c->fQ=0;
546
547
548 for (Int_t i=0; i<c->fMultiplicity; i++)
549 {
550 dig= (AliMUONdigit*)fDigits->UncheckedAt(c->fIndexMap[i]);
551 ix=dig->fPadX+c->fOffsetMap[i];
552 iy=dig->fPadY;
553 Int_t q=dig->fSignal;
554 if (dig->fPhysics >= dig->fSignal) {
555 c->fPhysicsMap[i]=2;
556 } else if (dig->fPhysics == 0) {
557 c->fPhysicsMap[i]=0;
558 } else c->fPhysicsMap[i]=1;
559//
560//
561// peak signal and track list
562 if (flag) {
563 if (q>c->fPeakSignal) {
564 c->fPeakSignal=q;
565/*
566 c->fTracks[0]=dig->fTracks[0];
567 c->fTracks[1]=dig->fTracks[1];
568 c->fTracks[2]=dig->fTracks[2];
569*/
570 //c->fTracks[0]=dig->fTrack;
571 c->fTracks[0]=dig->fHit;
572 c->fTracks[1]=dig->fTracks[0];
573 c->fTracks[2]=dig->fTracks[1];
574 }
575 } else {
576 if (c->fContMap[i] > frac) {
577 frac=c->fContMap[i];
578 c->fPeakSignal=q;
579/*
580 c->fTracks[0]=dig->fTracks[0];
581 c->fTracks[1]=dig->fTracks[1];
582 c->fTracks[2]=dig->fTracks[2];
583*/
584 //c->fTracks[0]=dig->fTrack;
585 c->fTracks[0]=dig->fHit;
586 c->fTracks[1]=dig->fTracks[0];
587 c->fTracks[2]=dig->fTracks[1];
588 }
589 }
590//
591 if (flag) {
592 fSegmentation->GetPadCxy(ix, iy, x, y);
593 c->fX += q*x;
594 c->fY += q*y;
595 c->fQ += q;
596 }
597
598 } // loop over digits
599
600 if (flag) {
601
602 c->fX/=c->fQ;
603 c->fX=fSegmentation->GetAnod(c->fX);
604 c->fY/=c->fQ;
605//
606// apply correction to the coordinate along the anode wire
607//
608 x=c->fX;
609 y=c->fY;
610 fSegmentation->GetPadIxy(x, y, ix, iy);
611 fSegmentation->GetPadCxy(ix, iy, x, y);
612 Int_t isec=fSegmentation->Sector(ix,iy);
613 TF1* CogCorr = fSegmentation->CorrFunc(isec-1);
614
615 if (CogCorr) {
616 Float_t YonPad=(c->fY-y)/fSegmentation->Dpy(isec);
617 c->fY=c->fY-CogCorr->Eval(YonPad, 0, 0);
618 }
619 }
620}
621
622
623void AliMUONClusterFinder::FindCluster(Int_t i, Int_t j, AliMUONRawCluster &c){
624//
625// Find clusters
626//
627//
628// Add i,j as element of the cluster
629//
630
631 Int_t idx = fHitMap->GetHitIndex(i,j);
632 AliMUONdigit* dig = (AliMUONdigit*) fHitMap->GetHit(i,j);
633 Int_t q=dig->fSignal;
634 if (q > TMath::Abs(c.fPeakSignal)) {
635 c.fPeakSignal=q;
636/*
637 c.fTracks[0]=dig->fTracks[0];
638 c.fTracks[1]=dig->fTracks[1];
639 c.fTracks[2]=dig->fTracks[2];
640*/
641 //c.fTracks[0]=dig->fTrack;
642 c.fTracks[0]=dig->fHit;
643 c.fTracks[1]=dig->fTracks[0];
644 c.fTracks[2]=dig->fTracks[1];
645 }
646//
647// Make sure that list of digits is ordered
648//
649 Int_t mu=c.fMultiplicity;
650 c.fIndexMap[mu]=idx;
651
652 if (dig->fPhysics >= dig->fSignal) {
653 c.fPhysicsMap[mu]=2;
654 } else if (dig->fPhysics == 0) {
655 c.fPhysicsMap[mu]=0;
656 } else c.fPhysicsMap[mu]=1;
657
658 if (mu > 0) {
659 for (Int_t ind=mu-1; ind>=0; ind--) {
660 Int_t ist=(c.fIndexMap)[ind];
661 Int_t ql=((AliMUONdigit*)fDigits
662 ->UncheckedAt(ist))->fSignal;
663 if (q>ql) {
664 c.fIndexMap[ind]=idx;
665 c.fIndexMap[ind+1]=ist;
666 } else {
667 break;
668 }
669 }
670 }
671
672 c.fMultiplicity++;
673
674 if (c.fMultiplicity >= 50 ) {
675 printf("FindCluster - multiplicity >50 %d \n",c.fMultiplicity);
676 c.fMultiplicity=49;
677 }
678
679// Prepare center of gravity calculation
680 Float_t x, y;
681 fSegmentation->GetPadCxy(i, j, x, y);
682 c.fX += q*x;
683 c.fY += q*y;
684 c.fQ += q;
685// Flag hit as taken
686 fHitMap->FlagHit(i,j);
687//
688// Now look recursively for all neighbours
689//
690 Int_t nn;
691 Int_t Xlist[kMaxNeighbours], Ylist[kMaxNeighbours];
692 fSegmentation->Neighbours(i,j,&nn,Xlist,Ylist);
693 for (Int_t in=0; in<nn; in++) {
694 Int_t ix=Xlist[in];
695 Int_t iy=Ylist[in];
696 if (fHitMap->TestHit(ix,iy)==unused) FindCluster(ix, iy, c);
697 }
698}
699
700//_____________________________________________________________________________
701
702void AliMUONClusterFinder::FindRawClusters()
703{
704 //
705 // simple MUON cluster finder from digits -- finds neighbours and
706 // fill the tree with raw clusters
707 //
e3a4d40e 708
a897a37a 709 if (!fNdigits) return;
710
711 fHitMap = new AliMUONHitMapA1(fSegmentation, fDigits);
712
e3a4d40e 713
a897a37a 714 AliMUONdigit *dig;
715
716 Int_t ndig;
717 Int_t nskip=0;
718 Int_t ncls=0;
719 fHitMap->FillHits();
720 for (ndig=0; ndig<fNdigits; ndig++) {
721 dig = (AliMUONdigit*)fDigits->UncheckedAt(ndig);
722 Int_t i=dig->fPadX;
723 Int_t j=dig->fPadY;
724 if (fHitMap->TestHit(i,j)==used ||fHitMap->TestHit(i,j)==empty) {
725 nskip++;
726 continue;
727 }
728 AliMUONRawCluster c;
729 c.fMultiplicity=0;
730 c.fPeakSignal=dig->fSignal;
731/*
732 c.fTracks[0]=dig->fTracks[0];
733 c.fTracks[1]=dig->fTracks[1];
734 c.fTracks[2]=dig->fTracks[2];
735*/
736 //c.fTracks[0]=dig->fTrack;
737 c.fTracks[0]=dig->fHit;
738 c.fTracks[1]=dig->fTracks[0];
739 c.fTracks[2]=dig->fTracks[1];
740 // tag the beginning of cluster list in a raw cluster
741 c.fNcluster[0]=-1;
742 FindCluster(i,j, c);
743 // center of gravity
744 c.fX /= c.fQ;
745 c.fX=fSegmentation->GetAnod(c.fX);
746 c.fY /= c.fQ;
747//
748// apply correction to the coordinate along the anode wire
749//
750 Int_t ix,iy;
751 Float_t x=c.fX;
752 Float_t y=c.fY;
753 fSegmentation->GetPadIxy(x, y, ix, iy);
754 fSegmentation->GetPadCxy(ix, iy, x, y);
755 Int_t isec=fSegmentation->Sector(ix,iy);
756 TF1* CogCorr=fSegmentation->CorrFunc(isec-1);
757 if (CogCorr) {
758 Float_t YonPad=(c.fY-y)/fSegmentation->Dpy(isec);
759 c.fY=c.fY-CogCorr->Eval(YonPad,0,0);
760 }
761
762//
763// Analyse cluster and decluster if necessary
764//
765 ncls++;
766 c.fNcluster[1]=fNRawClusters;
767 c.fClusterType=c.PhysicsContribution();
768 Decluster(&c);
769 fNPeaks=0;
770//
771//
772//
773// reset Cluster object
774 for (int k=0;k<c.fMultiplicity;k++) {
775 c.fIndexMap[k]=0;
776 }
777 c.fMultiplicity=0;
778 } // end loop ndig
779 delete fHitMap;
780}
781
782void AliMUONClusterFinder::
783CalibrateCOG()
784{
785 Float_t x[5];
786 Float_t y[5];
787 Int_t n, i;
788 TF1 func;
789 if (fSegmentation) {
790 fSegmentation->GiveTestPoints(n, x, y);
791 for (i=0; i<n; i++) {
792 Float_t xtest=x[i];
793 Float_t ytest=y[i];
794 SinoidalFit(xtest, ytest, func);
795 fSegmentation->SetCorrFunc(i, new TF1(func));
796 }
797 }
798}
799
800
801void AliMUONClusterFinder::
802SinoidalFit(Float_t x, Float_t y, TF1 &func)
803{
804//
805 static Int_t count=0;
806 char canvasname[3];
807 count++;
808 sprintf(canvasname,"c%d",count);
809
a6f39961 810 const Int_t ns=101;
a897a37a 811 Float_t xg[ns], yg[ns], xrg[ns], yrg[ns];
812 Float_t xsig[ns], ysig[ns];
813
814 AliMUONsegmentation *segmentation=fSegmentation;
815
816 Int_t ix,iy;
817 segmentation->GetPadIxy(x,y,ix,iy);
818 segmentation->GetPadCxy(ix,iy,x,y);
819 Int_t isec=segmentation->Sector(ix,iy);
820// Pad Limits
821 Float_t xmin = x-segmentation->Dpx(isec)/2;
822 Float_t ymin = y-segmentation->Dpy(isec)/2;
823//
824// Integration Limits
825 Float_t dxI=fResponse->SigmaIntegration()*fResponse->ChargeSpreadX();
826 Float_t dyI=fResponse->SigmaIntegration()*fResponse->ChargeSpreadY();
827
828//
829// Scanning
830//
831 Int_t i;
832 Float_t qp;
833//
834// y-position
835 Float_t yscan=ymin;
836 Float_t dy=segmentation->Dpy(isec)/(ns-1);
837
838 for (i=0; i<ns; i++) {
839//
840// Pad Loop
841//
842 Float_t sum=0;
843 Float_t qcheck=0;
844 segmentation->SigGenInit(x, yscan, 0);
845
846 for (segmentation->FirstPad(x, yscan, dxI, dyI);
847 segmentation->MorePads();
848 segmentation->NextPad())
849 {
850 qp=fResponse->IntXY(segmentation);
851 qp=TMath::Abs(qp);
852//
853//
854 if (qp > 1.e-4) {
855 qcheck+=qp;
856 Int_t ixs=segmentation->Ix();
857 Int_t iys=segmentation->Iy();
858 Float_t xs,ys;
859 segmentation->GetPadCxy(ixs,iys,xs,ys);
860 sum+=qp*ys;
861 }
862 } // Pad loop
863 Float_t ycog=sum/qcheck;
864 yg[i]=(yscan-y)/segmentation->Dpy(isec);
865 yrg[i]=(ycog-y)/segmentation->Dpy(isec);
866 ysig[i]=ycog-yscan;
867 yscan+=dy;
868 } // scan loop
869//
870// x-position
871 Float_t xscan=xmin;
872 Float_t dx=segmentation->Dpx(isec)/(ns-1);
873
874 for (i=0; i<ns; i++) {
875//
876// Pad Loop
877//
878 Float_t sum=0;
879 Float_t qcheck=0;
880 segmentation->SigGenInit(xscan, y, 0);
881
882 for (segmentation->FirstPad(xscan, y, dxI, dyI);
883 segmentation->MorePads();
884 segmentation->NextPad())
885 {
886 qp=fResponse->IntXY(segmentation);
887 qp=TMath::Abs(qp);
888//
889//
890 if (qp > 1.e-2) {
891 qcheck+=qp;
892 Int_t ixs=segmentation->Ix();
893 Int_t iys=segmentation->Iy();
894 Float_t xs,ys;
895 segmentation->GetPadCxy(ixs,iys,xs,ys);
896 sum+=qp*xs;
897 }
898 } // Pad loop
899 Float_t xcog=sum/qcheck;
900 xcog=segmentation->GetAnod(xcog);
901
902 xg[i]=(xscan-x)/segmentation->Dpx(isec);
903 xrg[i]=(xcog-x)/segmentation->Dpx(isec);
904 xsig[i]=xcog-xscan;
905 xscan+=dx;
906 }
907//
908// Creates a Root function based on function sinoid above
909// and perform the fit
910//
911 // TGraph *graphx = new TGraph(ns,xg ,xsig);
912 // TGraph *graphxr= new TGraph(ns,xrg,xsig);
913 // TGraph *graphy = new TGraph(ns,yg ,ysig);
914 TGraph *graphyr= new TGraph(ns,yrg,ysig);
915
916 Double_t sinoid(Double_t *x, Double_t *par);
917 new TF1("sinoidf",sinoid,0.5,0.5,5);
918 graphyr->Fit("sinoidf","Q");
919 func = *((TF1*)((graphyr->GetListOfFunctions())->At(0)));
920/*
921
922 TCanvas *c1=new TCanvas(canvasname,canvasname,400,10,600,700);
923 TPad* pad11 = new TPad("pad11"," ",0.01,0.51,0.49,0.99);
924 TPad* pad12 = new TPad("pad12"," ",0.51,0.51,0.99,0.99);
925 TPad* pad13 = new TPad("pad13"," ",0.01,0.01,0.49,0.49);
926 TPad* pad14 = new TPad("pad14"," ",0.51,0.01,0.99,0.49);
927 pad11->SetFillColor(11);
928 pad12->SetFillColor(11);
929 pad13->SetFillColor(11);
930 pad14->SetFillColor(11);
931 pad11->Draw();
932 pad12->Draw();
933 pad13->Draw();
934 pad14->Draw();
935
936//
937 pad11->cd();
938 graphx->SetFillColor(42);
939 graphx->SetMarkerColor(4);
940 graphx->SetMarkerStyle(21);
941 graphx->Draw("AC");
942 graphx->GetHistogram()->SetXTitle("x on pad");
943 graphx->GetHistogram()->SetYTitle("xcog-x");
944
945
946 pad12->cd();
947 graphxr->SetFillColor(42);
948 graphxr->SetMarkerColor(4);
949 graphxr->SetMarkerStyle(21);
950 graphxr->Draw("AP");
951 graphxr->GetHistogram()->SetXTitle("xcog on pad");
952 graphxr->GetHistogram()->SetYTitle("xcog-x");
953
954
955 pad13->cd();
956 graphy->SetFillColor(42);
957 graphy->SetMarkerColor(4);
958 graphy->SetMarkerStyle(21);
959 graphy->Draw("AF");
960 graphy->GetHistogram()->SetXTitle("y on pad");
961 graphy->GetHistogram()->SetYTitle("ycog-y");
962
963
964
965 pad14->cd();
966 graphyr->SetFillColor(42);
967 graphyr->SetMarkerColor(4);
968 graphyr->SetMarkerStyle(21);
969 graphyr->Draw("AF");
970 graphyr->GetHistogram()->SetXTitle("ycog on pad");
971 graphyr->GetHistogram()->SetYTitle("ycog-y");
972
973 c1->Update();
974*/
975}
976
977Double_t sinoid(Double_t *x, Double_t *par)
978{
979 Double_t arg = -2*TMath::Pi()*x[0];
980 Double_t fitval= par[0]*TMath::Sin(arg)+
981 par[1]*TMath::Sin(2*arg)+
982 par[2]*TMath::Sin(3*arg)+
983 par[3]*TMath::Sin(4*arg)+
984 par[4]*TMath::Sin(5*arg);
985 return fitval;
986 }
987
988
989Double_t DoubleGauss(Double_t *x, Double_t *par)
990{
991 Double_t arg1 = (x[0]-par[1])/0.18;
992 Double_t arg2 = (x[0]-par[3])/0.18;
993 Double_t fitval= par[0]*TMath::Exp(-arg1*arg1/2)
994 +par[2]*TMath::Exp(-arg2*arg2/2);
995 return fitval;
996 }
997
998Float_t DiscrCharge(Int_t i,Double_t *par)
999{
1000// par[0] x-position of first cluster
1001// par[1] y-position of first cluster
1002// par[2] x-position of second cluster
1003// par[3] y-position of second cluster
1004// par[4] charge fraction of first cluster
1005// 1-par[4] charge fraction of second cluster
1006
1007 static Float_t qtot;
1008 if (gFirst) {
1009 qtot=0;
1010 for (Int_t jbin=0; jbin<gNbins; jbin++) {
1011 qtot+=gCharge[jbin];
1012 }
1013 gFirst=0;
1014 //printf("\n sum of charge from DiscrCharge %f\n", qtot);
a6f39961 1015 gChargeTot=Int_t(qtot);
a897a37a 1016
1017 }
1018 gSegmentation->SetPad(gix[i], giy[i]);
1019// First Cluster
1020 gSegmentation->SetHit(par[0],par[1]);
1021 Float_t q1=gResponse->IntXY(gSegmentation);
1022
1023// Second Cluster
1024 gSegmentation->SetHit(par[2],par[3]);
1025 Float_t q2=gResponse->IntXY(gSegmentation);
1026
1027 Float_t value = qtot*(par[4]*q1+(1.-par[4])*q2);
1028 return value;
1029}
1030
1031//
1032// Minimisation function
e3a4d40e 1033void fcn(Int_t &, Double_t *, Double_t &f, Double_t *par, Int_t )
a897a37a 1034{
1035 Int_t i;
1036 Float_t delta;
1037 Float_t chisq=0;
1038 Float_t qcont=0;
1039 Float_t qtot=0;
1040
1041 for (i=0; i<gNbins; i++) {
1042 Float_t q0=gCharge[i];
1043 Float_t q1=DiscrCharge(i,par);
1044 delta=(q0-q1)/TMath::Sqrt(q0);
1045 chisq+=delta*delta;
1046 qcont+=q1;
1047 qtot+=q0;
1048 }
1049 chisq=chisq+=(qtot-qcont)*(qtot-qcont)*0.5;
1050 f=chisq;
1051}
1052
1053
1054
1055
1056
1057
1058
1059
1060
1061
1062