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