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