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