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