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