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