CRAB added
[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
8265fa96 16
17#include "AliRICHClusterFinder.h"
8265fa96 18#include "AliRun.h"
237c933d 19#include "AliRICH.h"
543d5224 20#include "AliRICHMap.h"
b251a2b5 21#include "AliRICHSDigit.h"
237c933d 22#include "AliRICHDigit.h"
23#include "AliRICHRawCluster.h"
543d5224 24#include "AliRICHParam.h"
237c933d 25#include <TTree.h>
8265fa96 26#include <TCanvas.h>
27#include <TH1.h>
488e98ba 28#include <TF1.h>
8265fa96 29#include <TPad.h>
30#include <TGraph.h>
bc808d18 31#include <TMinuit.h>
8265fa96 32
543d5224 33static AliSegmentation *gSegmentation;
8265fa96 34static AliRICHResponse* gResponse;
35static Int_t gix[500];
36static Int_t giy[500];
37static Float_t gCharge[500];
38static Int_t gNbins;
6ce834b4 39static Bool_t gFirst=kTRUE;
8265fa96 40static TMinuit *gMyMinuit ;
6ce834b4 41void fcn(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t);
8265fa96 42static Int_t gChargeTot;
43
44ClassImp(AliRICHClusterFinder)
543d5224 45//__________________________________________________________________________________________________
46AliRICHClusterFinder::AliRICHClusterFinder(AliRICH *pRICH)
47{//main ctor
48 Info("main ctor","Start.");
49
50 fRICH=pRICH;
51
52 fSegmentation=Rich()->C(1)->GetSegmentationModel();
53 fResponse =Rich()->C(1)->GetResponseModel();
8265fa96 54
543d5224 55 fDigits=0; fNdigits=0;
56 fChamber=0;
57 fHitMap=0;
58
59 fCogCorr = 0;
60 SetNperMax();
61 SetClusterSize();
543d5224 62 fNPeaks=-1;
63}//main ctor
64//__________________________________________________________________________________________________
8265fa96 65void AliRICHClusterFinder::Decluster(AliRICHRawCluster *cluster)
543d5224 66{// Decluster algorithm
67 Info("Decluster","Start.");
68 if(cluster->fMultiplicity==1||cluster->fMultiplicity==2){//Nothing special for 1- and 2-clusters
69 if(fNPeaks != 0) {cluster->fNcluster[0]=fNPeaks; cluster->fNcluster[1]=0;}
70 AddRawCluster(*cluster);
71 fNPeaks++;
72 }else if(cluster->fMultiplicity==3){// 3-cluster, check topology
543d5224 73 Centered(cluster);// ok, cluster is centered and added in Centered()
543d5224 74 }else{//4-and more-pad clusters
6ce834b4 75 if(cluster->fMultiplicity<= fMaxClusterSize){
543d5224 76 SplitByLocalMaxima(cluster);
543d5224 77 }//if <= fClusterSize
78 }//if multiplicity
79}//Decluster()
80//__________________________________________________________________________________________________
8265fa96 81Bool_t AliRICHClusterFinder::Centered(AliRICHRawCluster *cluster)
543d5224 82{//Is the cluster centered?
237c933d 83
543d5224 84 AliRICHDigit* dig;
85 dig= (AliRICHDigit*)fDigits->UncheckedAt(cluster->fIndexMap[0]);
86 Int_t x[kMaxNeighbours], y[kMaxNeighbours], xN[kMaxNeighbours], yN[kMaxNeighbours];
c712cb2f 87 Int_t nn=Rich()->Param()->PadNeighbours(dig->PadX(),dig->PadY(),x,y);
8265fa96 88
543d5224 89
90 Int_t nd=0;
91 for (Int_t i=0; i<nn; i++){//neighbours loop
92 if(fHitMap->TestHit(x[i],y[i]) == kUsed){
93 xN[nd]=x[i];
94 yN[nd]=y[i];
95 nd++;
8265fa96 96 }
543d5224 97 }//neighbours loop
98
99 if(nd==2){// cluster is centered !
8265fa96 100 if (fNPeaks != 0) {
101 cluster->fNcluster[0]=fNPeaks;
102 cluster->fNcluster[1]=0;
103 }
104 cluster->fCtype=0;
105 AddRawCluster(*cluster);
106 fNPeaks++;
107 return kTRUE;
108 } else if (nd ==1) {
8265fa96 109// Highest signal on an edge, split cluster into 2+1
543d5224 110// who is the neighbour ?
237c933d 111 Int_t nind=fHitMap->GetHitIndex(xN[0], yN[0]);
8265fa96 112 Int_t i1= (nind==cluster->fIndexMap[1]) ? 1:2;
113 Int_t i2= (nind==cluster->fIndexMap[1]) ? 2:1;
8265fa96 114// 2-cluster
115 AliRICHRawCluster cnew;
116 if (fNPeaks == 0) {
117 cnew.fNcluster[0]=-1;
118 cnew.fNcluster[1]=fNRawClusters;
119 } else {
120 cnew.fNcluster[0]=fNPeaks;
121 cnew.fNcluster[1]=0;
122 }
123 cnew.fMultiplicity=2;
124 cnew.fIndexMap[0]=cluster->fIndexMap[0];
125 cnew.fIndexMap[1]=cluster->fIndexMap[i1];
126 FillCluster(&cnew);
127 cnew.fClusterType=cnew.PhysicsContribution();
128 AddRawCluster(cnew);
129 fNPeaks++;
8265fa96 130// 1-cluster
131 cluster->fMultiplicity=1;
132 cluster->fIndexMap[0]=cluster->fIndexMap[i2];
133 cluster->fIndexMap[1]=0;
134 cluster->fIndexMap[2]=0;
135 FillCluster(cluster);
136 if (fNPeaks != 0) {
137 cluster->fNcluster[0]=fNPeaks;
138 cluster->fNcluster[1]=0;
139 }
140 cluster->fClusterType=cluster->PhysicsContribution();
141 AddRawCluster(*cluster);
142 fNPeaks++;
143 return kFALSE;
144 } else {
543d5224 145 Warning("Centered","\n Completely screwed up %d !! \n",nd);
8265fa96 146
543d5224 147 }
148 return kFALSE;
149}//Centered()
150//__________________________________________________________________________________________________
8265fa96 151void AliRICHClusterFinder::SplitByLocalMaxima(AliRICHRawCluster *c)
543d5224 152{// Split the cluster according to the number of maxima inside
6ce834b4 153 Info("SplitbyLocalMaxima","Start.");
154
155 AliRICHDigit* dig[100], *digt;
8265fa96 156 Int_t ix[100], iy[100], q[100];
c712cb2f 157 Double_t x[100], y[100];
8265fa96 158 Int_t i; // loops over digits
159 Int_t j; // loops over local maxima
8265fa96 160 Int_t mul=c->fMultiplicity;
8265fa96 161// dump digit information into arrays
543d5224 162 for (i=0; i<mul; i++){
163 dig[i]= (AliRICHDigit*)fDigits->UncheckedAt(c->fIndexMap[i]);
164 ix[i]= dig[i]->PadX();
165 iy[i]= dig[i]->PadY();
166 q[i] = dig[i]->Signal();
c712cb2f 167 AliRICHParam::Pad2Loc(ix[i], iy[i], x[i], y[i]);
543d5224 168 }
8265fa96 169// Find local maxima
237c933d 170 Bool_t isLocal[100];
171 Int_t nLocal=0;
172 Int_t associatePeak[100];
173 Int_t indLocal[100];
8265fa96 174 Int_t nn;
237c933d 175 Int_t xNei[kMaxNeighbours], yNei[kMaxNeighbours];
8265fa96 176 for (i=0; i<mul; i++) {
237c933d 177 fSegmentation->Neighbours(ix[i], iy[i], &nn, xNei, yNei);
178 isLocal[i]=kTRUE;
8265fa96 179 for (j=0; j<nn; j++) {
237c933d 180 if (fHitMap->TestHit(xNei[j], yNei[j])==kEmpty) continue;
181 digt=(AliRICHDigit*) fHitMap->GetHit(xNei[j], yNei[j]);
6e585aa2 182 if (digt->Signal() > q[i]) {
237c933d 183 isLocal[i]=kFALSE;
8265fa96 184 break;
8265fa96 185// handle special case of neighbouring pads with equal signal
6e585aa2 186 } else if (digt->Signal() == q[i]) {
237c933d 187 if (nLocal >0) {
188 for (Int_t k=0; k<nLocal; k++) {
189 if (xNei[j]==ix[indLocal[k]] && yNei[j]==iy[indLocal[k]]){
190 isLocal[i]=kFALSE;
8265fa96 191 }
192 }
193 }
194 }
195 } // loop over next neighbours
196 // Maxima should not be on the edge
237c933d 197 if (isLocal[i]) {
198 indLocal[nLocal]=i;
199 nLocal++;
8265fa96 200 }
201 } // loop over all digits
543d5224 202// If only one local maximum found but multiplicity is high take global maximum from the list of digits.
237c933d 203 if (nLocal==1 && mul>5) {
8265fa96 204 Int_t nnew=0;
205 for (i=0; i<mul; i++) {
237c933d 206 if (!isLocal[i]) {
207 indLocal[nLocal]=i;
208 isLocal[i]=kTRUE;
209 nLocal++;
8265fa96 210 nnew++;
211 }
212 if (nnew==1) break;
213 }
214 }
6ce834b4 215 if(nLocal==2) {// If number of local maxima is 2 try to fit a double gaussian
216
8265fa96 217// Initialise global variables for fit
6ce834b4 218 gFirst=kTRUE;
8265fa96 219 gSegmentation=fSegmentation;
220 gResponse =fResponse;
221 gNbins=mul;
222
223 for (i=0; i<mul; i++) {
224 gix[i]=ix[i];
225 giy[i]=iy[i];
226 gCharge[i]=Float_t(q[i]);
227 }
6ce834b4 228 if (gFirst) gMyMinuit = new TMinuit(5);
229
8265fa96 230 gMyMinuit->SetFCN(fcn);
231 gMyMinuit->mninit(5,10,7);
232 Double_t arglist[20];
8265fa96 233 arglist[0]=1;
8265fa96 234// Set starting values
235 static Double_t vstart[5];
237c933d 236 vstart[0]=x[indLocal[0]];
237 vstart[1]=y[indLocal[0]];
238 vstart[2]=x[indLocal[1]];
239 vstart[3]=y[indLocal[1]];
6ce834b4 240 vstart[4]=Float_t(q[indLocal[0]])/Float_t(q[indLocal[0]]+q[indLocal[1]]);
8265fa96 241// lower and upper limits
242 static Double_t lower[5], upper[5];
6ce834b4 243 lower[0]=vstart[0]-AliRICHParam::PadSizeX()/2;
244 lower[1]=vstart[1]-AliRICHParam::PadSizeY()/2;
8265fa96 245
6ce834b4 246 upper[0]=vstart[0]+AliRICHParam::PadSizeX()/2;
247 upper[1]=vstart[1]+AliRICHParam::PadSizeY()/2;
8265fa96 248
6ce834b4 249 lower[2]=vstart[2]-AliRICHParam::PadSizeX()/2;
250 lower[3]=vstart[3]-AliRICHParam::PadSizeY()/2;
8265fa96 251
6ce834b4 252 upper[2]=vstart[2]+AliRICHParam::PadSizeX()/2;
253 upper[3]=vstart[3]+AliRICHParam::PadSizeY()/2;
8265fa96 254
255 lower[4]=0.;
256 upper[4]=1.;
257// step sizes
258 static Double_t step[5]={0.005, 0.03, 0.005, 0.03, 0.01};
6ce834b4 259 Int_t iErr;
8265fa96 260
6ce834b4 261 gMyMinuit->mnparm(0,"x1",vstart[0],step[0],lower[0],upper[0],iErr);
262 gMyMinuit->mnparm(1,"y1",vstart[1],step[1],lower[1],upper[1],iErr);
263 gMyMinuit->mnparm(2,"x2",vstart[2],step[2],lower[2],upper[2],iErr);
264 gMyMinuit->mnparm(3,"y2",vstart[3],step[3],lower[3],upper[3],iErr);
265 gMyMinuit->mnparm(4,"a0",vstart[4],step[4],lower[4],upper[4],iErr);
8265fa96 266// ready for minimisation
267 gMyMinuit->SetPrintLevel(-1);
6ce834b4 268 gMyMinuit->mnexcm("SET OUT", arglist, 0, iErr);
8265fa96 269 arglist[0]= -1;
270 arglist[1]= 0;
271
6ce834b4 272 gMyMinuit->mnexcm("SET NOGR", arglist, 0, iErr);
273 gMyMinuit->mnexcm("SIMPLEX", arglist, 0, iErr);
274 gMyMinuit->mnexcm("MIGRAD", arglist, 0, iErr);
275 gMyMinuit->mnexcm("EXIT" , arglist, 0, iErr);
8265fa96 276
277 Double_t xrec[2], yrec[2], qfrac;
278 TString chname;
279 Double_t epxz, b1, b2;
6ce834b4 280 gMyMinuit->mnpout(0, chname, xrec[0], epxz, b1, b2, iErr);
281 gMyMinuit->mnpout(1, chname, yrec[0], epxz, b1, b2, iErr);
282 gMyMinuit->mnpout(2, chname, xrec[1], epxz, b1, b2, iErr);
283 gMyMinuit->mnpout(3, chname, yrec[1], epxz, b1, b2, iErr);
284 gMyMinuit->mnpout(4, chname, qfrac, epxz, b1, b2, iErr);
285
286 cout<<"xrex[0]="<<xrec[0]<<"yrec[0]="<<yrec[0]<<"xrec[1]="<<xrec[1]<<"yrec[1]="<<yrec[1]<<"qfrac="<<qfrac<<endl;
287 for (j=0; j<2; j++) { // One cluster for each maximum
8265fa96 288 AliRICHRawCluster cnew;
289 if (fNPeaks == 0) {
290 cnew.fNcluster[0]=-1;
291 cnew.fNcluster[1]=fNRawClusters;
292 } else {
293 cnew.fNcluster[0]=fNPeaks;
294 cnew.fNcluster[1]=0;
295 }
296 cnew.fMultiplicity=0;
297 cnew.fX=Float_t(xrec[j]);
298 cnew.fY=Float_t(yrec[j]);
299 if (j==0) {
300 cnew.fQ=Int_t(gChargeTot*qfrac);
301 } else {
302 cnew.fQ=Int_t(gChargeTot*(1-qfrac));
303 }
8265fa96 304 for (i=0; i<mul; i++) {
305 cnew.fIndexMap[cnew.fMultiplicity]=c->fIndexMap[i];
c712cb2f 306 TVector3 x3(xrec[j],yrec[j],0);
307 cnew.fContMap[cnew.fMultiplicity]=AliRICHParam::Loc2PadFrac(x3,gix[i], giy[i]);
8265fa96 308 cnew.fMultiplicity++;
309 }
310 FillCluster(&cnew,0);
8265fa96 311 cnew.fClusterType=cnew.PhysicsContribution();
312 AddRawCluster(cnew);
313 fNPeaks++;
314 }
6ce834b4 315 }//if 2 maximum in cluster
8265fa96 316 Bool_t fitted=kTRUE;
317
6ce834b4 318 if (nLocal >2 || !fitted) {
543d5224 319 // Check if enough local clusters have been found, if not add global maxima to the list
8265fa96 320 Int_t nPerMax;
237c933d 321 if (nLocal!=0) {
322 nPerMax=mul/nLocal;
8265fa96 323 } else {
543d5224 324 Warning("SplitByLocalMaxima","no local maximum found");
8265fa96 325 nPerMax=fNperMax+1;
326 }
327
328 if (nPerMax > fNperMax) {
237c933d 329 Int_t nGlob=mul/fNperMax-nLocal+1;
8265fa96 330 if (nGlob > 0) {
331 Int_t nnew=0;
332 for (i=0; i<mul; i++) {
237c933d 333 if (!isLocal[i]) {
334 indLocal[nLocal]=i;
335 isLocal[i]=kTRUE;
336 nLocal++;
8265fa96 337 nnew++;
338 }
339 if (nnew==nGlob) break;
340 }
341 }
342 }
543d5224 343 for (i=0; i<mul; i++) { // Associate hits to peaks
8265fa96 344 Float_t dmin=1.E10;
345 Float_t qmax=0;
237c933d 346 if (isLocal[i]) continue;
347 for (j=0; j<nLocal; j++) {
348 Int_t il=indLocal[j];
8265fa96 349 Float_t d=TMath::Sqrt((x[i]-x[il])*(x[i]-x[il])
350 +(y[i]-y[il])*(y[i]-y[il]));
351 Float_t ql=q[il];
543d5224 352 if (d<dmin) { // Select nearest peak
8265fa96 353 dmin=d;
354 qmax=ql;
237c933d 355 associatePeak[i]=j;
543d5224 356 } else if (d==dmin) { // If more than one take highest peak
8265fa96 357 if (ql>qmax) {
358 dmin=d;
359 qmax=ql;
237c933d 360 associatePeak[i]=j;
8265fa96 361 }
362 }
363 }
543d5224 364 }
8265fa96 365 // One cluster for each maximum
237c933d 366 for (j=0; j<nLocal; j++) {
8265fa96 367 AliRICHRawCluster cnew;
368 if (fNPeaks == 0) {
369 cnew.fNcluster[0]=-1;
370 cnew.fNcluster[1]=fNRawClusters;
371 } else {
372 cnew.fNcluster[0]=fNPeaks;
373 cnew.fNcluster[1]=0;
374 }
237c933d 375 cnew.fIndexMap[0]=c->fIndexMap[indLocal[j]];
8265fa96 376 cnew.fMultiplicity=1;
377 for (i=0; i<mul; i++) {
237c933d 378 if (isLocal[i]) continue;
379 if (associatePeak[i]==j) {
8265fa96 380 cnew.fIndexMap[cnew.fMultiplicity]=c->fIndexMap[i];
381 cnew.fMultiplicity++;
382 }
383 }
384 FillCluster(&cnew);
385 cnew.fClusterType=cnew.PhysicsContribution();
386 AddRawCluster(cnew);
387 fNPeaks++;
388 }
389 }
543d5224 390}//SplitByLocalMaxima(AliRICHRawCluster *c)
391//__________________________________________________________________________________________________
8265fa96 392void AliRICHClusterFinder::FillCluster(AliRICHRawCluster* c, Int_t flag)
543d5224 393{// Completes cluster information starting from list of digits
8265fa96 394 AliRICHDigit* dig;
c712cb2f 395 Double_t x, y;
8265fa96 396 Int_t ix, iy;
6ce834b4 397 Float_t fraction=0;
8265fa96 398
399 c->fPeakSignal=0;
400 if (flag) {
401 c->fX=0;
402 c->fY=0;
403 c->fQ=0;
404 }
8265fa96 405
406
543d5224 407 for (Int_t i=0; i<c->fMultiplicity; i++){
8265fa96 408 dig= (AliRICHDigit*)fDigits->UncheckedAt(c->fIndexMap[i]);
6ce834b4 409 ix=dig->PadX();
6e585aa2 410 iy=dig->PadY();
411 Int_t q=dig->Signal();
412 if (dig->Physics() >= dig->Signal()) {
8265fa96 413 c->fPhysicsMap[i]=2;
6e585aa2 414 } else if (dig->Physics() == 0) {
8265fa96 415 c->fPhysicsMap[i]=0;
416 } else c->fPhysicsMap[i]=1;
8265fa96 417// peak signal and track list
418 if (flag) {
419 if (q>c->fPeakSignal) {
420 c->fPeakSignal=q;
6e585aa2 421 c->fTracks[0]=dig->Hit();
422 c->fTracks[1]=dig->Track(0);
423 c->fTracks[2]=dig->Track(1);
8265fa96 424 }
425 } else {
6ce834b4 426 if (c->fContMap[i] > fraction) {
427 fraction=c->fContMap[i];
8265fa96 428 c->fPeakSignal=q;
6e585aa2 429 c->fTracks[0]=dig->Hit();
430 c->fTracks[1]=dig->Track(0);
431 c->fTracks[2]=dig->Track(1);
8265fa96 432 }
433 }
8265fa96 434 if (flag) {
c712cb2f 435 AliRICHParam::Pad2Loc(ix,iy,x,y);
8265fa96 436 c->fX += q*x;
437 c->fY += q*y;
438 c->fQ += q;
439 }
440
441 } // loop over digits
442
443 if (flag) {
444
445 c->fX/=c->fQ;
446 c->fX=fSegmentation->GetAnod(c->fX);
447 c->fY/=c->fQ;
8265fa96 448// apply correction to the coordinate along the anode wire
8265fa96 449 x=c->fX;
450 y=c->fY;
c712cb2f 451 AliRICHParam::Loc2Pad(x,y,ix,iy);
452 AliRICHParam::Pad2Loc(ix,iy,x,y);
8265fa96 453 Int_t isec=fSegmentation->Sector(ix,iy);
237c933d 454 TF1* cogCorr = fSegmentation->CorrFunc(isec-1);
8265fa96 455
237c933d 456 if (cogCorr) {
457 Float_t yOnPad=(c->fY-y)/fSegmentation->Dpy(isec);
458 c->fY=c->fY-cogCorr->Eval(yOnPad, 0, 0);
8265fa96 459 }
460 }
6ce834b4 461}//FillCluster()
543d5224 462//__________________________________________________________________________________________________
463void AliRICHClusterFinder::AddDigit2Cluster(Int_t i, Int_t j, AliRICHRawCluster &c)
464{//Find clusters Add i,j as element of the cluster
465 Info("AddDigit2Cluster","Start with digit(%i,%i)",i,j);
466
467 Int_t idx = fHitMap->GetHitIndex(i,j);
468 AliRICHDigit* dig = (AliRICHDigit*) fHitMap->GetHit(i,j);
469 Int_t q=dig->Signal();
470 if(q>TMath::Abs(c.fPeakSignal)){
8265fa96 471 c.fPeakSignal=q;
6e585aa2 472 c.fTracks[0]=dig->Hit();
473 c.fTracks[1]=dig->Track(0);
474 c.fTracks[2]=dig->Track(1);
8265fa96 475 }
8265fa96 476// Make sure that list of digits is ordered
8265fa96 477 Int_t mu=c.fMultiplicity;
478 c.fIndexMap[mu]=idx;
479
6e585aa2 480 if (dig->Physics() >= dig->Signal()) {
8265fa96 481 c.fPhysicsMap[mu]=2;
6e585aa2 482 } else if (dig->Physics() == 0) {
8265fa96 483 c.fPhysicsMap[mu]=0;
484 } else c.fPhysicsMap[mu]=1;
485
486 if (mu > 0) {
487 for (Int_t ind=mu-1; ind>=0; ind--) {
488 Int_t ist=(c.fIndexMap)[ind];
543d5224 489 Int_t ql=((AliRICHDigit*)fDigits->UncheckedAt(ist))->Signal();
8265fa96 490 if (q>ql) {
491 c.fIndexMap[ind]=idx;
492 c.fIndexMap[ind+1]=ist;
493 } else {
494 break;
495 }
496 }
497 }
498
543d5224 499 c.fMultiplicity++;
8265fa96 500 if (c.fMultiplicity >= 50 ) {
543d5224 501 Info("AddDigit2CLuster","multiplicity >50 %d \n",c.fMultiplicity);
8265fa96 502 c.fMultiplicity=49;
503 }
c712cb2f 504 Double_t x,y;// Prepare center of gravity calculation
505 AliRICHParam::Pad2Loc(i,j,x,y);
543d5224 506 c.fX+=q*x; c.fY+=q*y; c.fQ += q;
507 fHitMap->FlagHit(i,j);// Flag hit as taken
8265fa96 508
8265fa96 509
543d5224 510 Int_t xList[4], yList[4]; // Now look recursively for all neighbours
c712cb2f 511 for (Int_t iNei=0;iNei<Rich()->Param()->PadNeighbours(i,j,xList,yList);iNei++)
543d5224 512 if(fHitMap->TestHit(xList[iNei],yList[iNei])==kUnused) AddDigit2Cluster(xList[iNei],yList[iNei],c);
513}//AddDigit2Cluster()
514//__________________________________________________________________________________________________
8265fa96 515void AliRICHClusterFinder::FindRawClusters()
543d5224 516{//finds neighbours and fill the tree with raw clusters
517 Info("FindRawClusters","Start for Chamber %i.",fChamber);
518
519 if(!fNdigits)return;
8265fa96 520
543d5224 521 fHitMap=new AliRICHMap(fDigits);
8265fa96 522
543d5224 523 for(Int_t iDigN=0;iDigN<fNdigits;iDigN++){//digits loop
524 AliRICHDigit *dig=(AliRICHDigit*)fDigits->UncheckedAt(iDigN);
525 Int_t i=dig->PadX(); Int_t j=dig->PadY();
526 if(fHitMap->TestHit(i,j)==kUsed||fHitMap->TestHit(i,j)==kEmpty) continue;
a2f7eaf6 527
543d5224 528 AliRICHRawCluster c;
529 c.fMultiplicity=0; c.fPeakSignal=dig->Signal();
530 c.fTracks[0]=dig->Hit();c.fTracks[1]=dig->Track(0);c.fTracks[2]=dig->Track(1);
531 c.fNcluster[0]=-1;// tag the beginning of cluster list in a raw cluster
532
533 AddDigit2Cluster(i,j,c);//form initial cluster
534
535 c.fX /= c.fQ; // center of gravity
536 //c.fX=fSegmentation->GetAnod(c.fX);
537 c.fY /= c.fQ;
6ce834b4 538 //AddRawCluster(c);
543d5224 539
540// Int_t ix,iy;// apply correction to the coordinate along the anode wire
541// Float_t x=c.fX, y=c.fY;
c712cb2f 542// Rich()->Param()->Loc2Pad(x,y,ix,iy);
543// Rich()->Param()->Pad2Loc(ix,iy,x,y);
543d5224 544// Int_t isec=fSegmentation->Sector(ix,iy);
545// TF1* cogCorr=fSegmentation->CorrFunc(isec-1);
546// if(cogCorr){
547// Float_t yOnPad=(c.fY-y)/fSegmentation->Dpy(isec);
548// c.fY=c.fY-cogCorr->Eval(yOnPad,0,0);
549// }
550
551 c.fNcluster[1]=fNRawClusters; c.fClusterType=c.PhysicsContribution();
552
8265fa96 553 Decluster(&c);
543d5224 554
8265fa96 555 fNPeaks=0;
8265fa96 556
543d5224 557 c.fMultiplicity=0; for(int k=0;k<c.fMultiplicity;k++) c.fIndexMap[k]=0;//reset cluster object
558 }//digits loop
559 delete fHitMap;
560 Info("FindRawClusters","Stop.");
561}//FindRawClusters()
562//__________________________________________________________________________________________________
563void AliRICHClusterFinder::CalibrateCOG()
564{// Calibration
237c933d 565
8265fa96 566 Float_t x[5];
567 Float_t y[5];
568 Int_t n, i;
8265fa96 569 if (fSegmentation) {
bc808d18 570 TF1 *func;
8265fa96 571 fSegmentation->GiveTestPoints(n, x, y);
572 for (i=0; i<n; i++) {
bc808d18 573 func = 0;
8265fa96 574 Float_t xtest=x[i];
575 Float_t ytest=y[i];
576 SinoidalFit(xtest, ytest, func);
bc808d18 577 if (func) fSegmentation->SetCorrFunc(i, new TF1(*func));
8265fa96 578 }
579 }
543d5224 580}//CalibrateCOG()
581//__________________________________________________________________________________________________
c712cb2f 582void AliRICHClusterFinder::SinoidalFit(Double_t x, Double_t y, TF1 *func)
543d5224 583{//Sinoidal fit
584 static Int_t count=0;
a2f7eaf6 585
8265fa96 586 count++;
8265fa96 587
237c933d 588 const Int_t kNs=101;
589 Float_t xg[kNs], yg[kNs], xrg[kNs], yrg[kNs];
590 Float_t xsig[kNs], ysig[kNs];
8265fa96 591
8265fa96 592 Int_t ix,iy;
c712cb2f 593 AliRICHParam::Loc2Pad(x,y,ix,iy);
594 AliRICHParam::Pad2Loc(ix,iy,x,y);
543d5224 595 Int_t isec=fSegmentation->Sector(ix,iy);
8265fa96 596// Pad Limits
6ce834b4 597 Float_t xmin = x-Rich()->Param()->PadSizeX()/2;
598 Float_t ymin = y-Rich()->Param()->PadSizeY()/2;
8265fa96 599//
c712cb2f 600// Integration Limits
601 Float_t dxI=Rich()->Param()->MathiensonDeltaX();
602 Float_t dyI=Rich()->Param()->MathiensonDeltaY();
8265fa96 603
604//
605// Scanning
606//
607 Int_t i;
543d5224 608 Float_t qp=0;
609
8265fa96 610// y-position
611 Float_t yscan=ymin;
6ce834b4 612 Float_t dy=Rich()->Param()->PadSizeY()/(kNs-1);
8265fa96 613
543d5224 614 for (i=0; i<kNs; i++) {// Pad Loop
8265fa96 615 Float_t sum=0;
616 Float_t qcheck=0;
543d5224 617 fSegmentation->SigGenInit(x, yscan, 0);
8265fa96 618
543d5224 619 for (fSegmentation->FirstPad(x, yscan,0, dxI, dyI);
620 fSegmentation->MorePads();
621 fSegmentation->NextPad())
8265fa96 622 {
543d5224 623 qp=fResponse->IntXY(fSegmentation);
8265fa96 624 qp=TMath::Abs(qp);
8265fa96 625 if (qp > 1.e-4) {
626 qcheck+=qp;
543d5224 627 Int_t ixs=fSegmentation->Ix();
628 Int_t iys=fSegmentation->Iy();
c712cb2f 629 Double_t xs,ys;
630 AliRICHParam::Pad2Loc(ixs,iys,xs,ys);
8265fa96 631 sum+=qp*ys;
632 }
633 } // Pad loop
634 Float_t ycog=sum/qcheck;
543d5224 635 yg[i]=(yscan-y)/fSegmentation->Dpy(isec);
636 yrg[i]=(ycog-y)/fSegmentation->Dpy(isec);
8265fa96 637 ysig[i]=ycog-yscan;
638 yscan+=dy;
639 } // scan loop
8265fa96 640// x-position
641 Float_t xscan=xmin;
543d5224 642 Float_t dx=fSegmentation->Dpx(isec)/(kNs-1);
8265fa96 643
543d5224 644 for (i=0; i<kNs; i++) {// Pad Loop
8265fa96 645 Float_t sum=0;
646 Float_t qcheck=0;
543d5224 647 fSegmentation->SigGenInit(xscan, y, 0);
8265fa96 648
543d5224 649 for (fSegmentation->FirstPad(xscan, y, 0, dxI, dyI);
650 fSegmentation->MorePads();
651 fSegmentation->NextPad())
8265fa96 652 {
543d5224 653 qp=fResponse->IntXY(fSegmentation);
8265fa96 654 qp=TMath::Abs(qp);
8265fa96 655 if (qp > 1.e-2) {
656 qcheck+=qp;
543d5224 657 Int_t ixs=fSegmentation->Ix();
658 Int_t iys=fSegmentation->Iy();
c712cb2f 659 Double_t xs,ys;
660 AliRICHParam::Pad2Loc(ixs,iys,xs,ys);
8265fa96 661 sum+=qp*xs;
662 }
663 } // Pad loop
664 Float_t xcog=sum/qcheck;
543d5224 665 xcog=fSegmentation->GetAnod(xcog);
8265fa96 666
543d5224 667 xg[i]=(xscan-x)/fSegmentation->Dpx(isec);
668 xrg[i]=(xcog-x)/fSegmentation->Dpx(isec);
8265fa96 669 xsig[i]=xcog-xscan;
670 xscan+=dx;
671 }
543d5224 672// Creates a Root function based on function sinoid above and perform the fit
237c933d 673 TGraph *graphyr= new TGraph(kNs,yrg,ysig);
8265fa96 674 Double_t sinoid(Double_t *x, Double_t *par);
675 new TF1("sinoidf",sinoid,0.5,0.5,5);
676 graphyr->Fit("sinoidf","Q");
bc808d18 677 func = (TF1*)graphyr->GetListOfFunctions()->At(0);
543d5224 678}//SinoidalFit()
679//__________________________________________________________________________________________________
8265fa96 680Double_t sinoid(Double_t *x, Double_t *par)
543d5224 681{// Sinoid function
237c933d 682
8265fa96 683 Double_t arg = -2*TMath::Pi()*x[0];
684 Double_t fitval= par[0]*TMath::Sin(arg)+
685 par[1]*TMath::Sin(2*arg)+
686 par[2]*TMath::Sin(3*arg)+
687 par[3]*TMath::Sin(4*arg)+
688 par[4]*TMath::Sin(5*arg);
689 return fitval;
543d5224 690}//sinoid()
691//__________________________________________________________________________________________________
8265fa96 692Double_t DoubleGauss(Double_t *x, Double_t *par)
543d5224 693{//Double gaussian function
694 Double_t arg1 = (x[0]-par[1])/0.18;
695 Double_t arg2 = (x[0]-par[3])/0.18;
696 return par[0]*TMath::Exp(-arg1*arg1/2)+par[2]*TMath::Exp(-arg2*arg2/2);
697}
698//__________________________________________________________________________________________________
8265fa96 699Float_t DiscrCharge(Int_t i,Double_t *par)
700{
701// par[0] x-position of first cluster
702// par[1] y-position of first cluster
703// par[2] x-position of second cluster
704// par[3] y-position of second cluster
705// par[4] charge fraction of first cluster
706// 1-par[4] charge fraction of second cluster
707
708 static Float_t qtot;
709 if (gFirst) {
710 qtot=0;
711 for (Int_t jbin=0; jbin<gNbins; jbin++) {
712 qtot+=gCharge[jbin];
713 }
6ce834b4 714 gFirst=kFALSE;
8265fa96 715 gChargeTot=Int_t(qtot);
716
717 }
c712cb2f 718 TVector3 x3(par[0],par[1],0);
719 Float_t q1=AliRICHParam::Loc2PadFrac(x3,gix[i],giy[i]);
720 x3.SetX(par[2]);x3.SetY(par[3]);
721 Float_t q2=AliRICHParam::Loc2PadFrac(x3,gix[i],giy[i]);
6ce834b4 722// cout<<"qtot="<<gChargeTot<<" q1="<<q1<<" q2="<<q2<<" px="<<gix[i]<<" py="<<giy[i]<<endl;
8265fa96 723
724 Float_t value = qtot*(par[4]*q1+(1.-par[4])*q2);
725 return value;
543d5224 726}//DiscrCharge(Int_t i,Double_t *par)
727//__________________________________________________________________________________________________
6ce834b4 728void fcn(Int_t &npar, Double_t */*gin*/, Double_t &f, Double_t *par, Int_t)
543d5224 729{// Minimisation function
cc23c5c6 730 npar=1;
8265fa96 731 Int_t i;
732 Float_t delta;
733 Float_t chisq=0;
734 Float_t qcont=0;
735 Float_t qtot=0;
736
737 for (i=0; i<gNbins; i++) {
738 Float_t q0=gCharge[i];
739 Float_t q1=DiscrCharge(i,par);
740 delta=(q0-q1)/TMath::Sqrt(q0);
741 chisq+=delta*delta;
742 qcont+=q1;
743 qtot+=q0;
744 }
6ce834b4 745// chisq=chisq+=(qtot-qcont)*(qtot-qcont)*0.5;
8265fa96 746 f=chisq;
543d5224 747}//
748//__________________________________________________________________________________________________
749void AliRICHClusterFinder::Exec()
237c933d 750{
543d5224 751 Info("Exec","Start.");
752
753 Rich()->GetLoader()->LoadDigits();
754
755 for(Int_t iEventN=0;iEventN<gAlice->GetEventsPerRun();iEventN++){//events loop
756 gAlice->GetRunLoader()->GetEvent(iEventN);
237c933d 757
543d5224 758 Rich()->GetLoader()->MakeTree("R"); Rich()->MakeBranch("R");
759 Rich()->ResetDigitsOld(); Rich()->ResetRawClusters();
760
761 Rich()->GetLoader()->TreeD()->GetEntry(0);
762 for(fChamber=1;fChamber<=kNCH;fChamber++){//chambers loop
763 fDigits=Rich()->DigitsOld(fChamber); fNdigits=fDigits->GetEntries();
764
765 FindRawClusters();
766
767 }//chambers loop
768
769 Rich()->GetLoader()->TreeR()->Fill();
770 Rich()->GetLoader()->WriteRecPoints("OVERWRITE");
771 }//events loop
772 Rich()->GetLoader()->UnloadDigits(); Rich()->GetLoader()->UnloadRecPoints();
773 Rich()->ResetDigitsOld(); Rich()->ResetRawClusters();
774 Info("Exec","Stop.");
775}//Exec()
776//__________________________________________________________________________________________________