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