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