]> git.uio.no Git - u/mrichter/AliRoot.git/blame - ITS/AliITSclustererV2.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
New member function RecPoints2Clusters
[u/mrichter/AliRoot.git] / ITS / AliITSclustererV2.cxx
CommitLineData
7f4044f0 1//-------------------------------------------------------------------------
2// Implementation of the ITS clusterer V2 class
3//
4// Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch
5//-------------------------------------------------------------------------
6//uncomment the line below for running with V1 cluster finder classes
7//#define V1
8
9#include "AliRun.h"
10
11#include "AliITSclustererV2.h"
12#include "AliITSclusterV2.h"
13
14#include <Riostream.h>
15#include <TFile.h>
16#include <TTree.h>
17#include <TClonesArray.h>
18#include "AliITSgeom.h"
19#include "AliITSdigit.h"
20
21ClassImp(AliITSclustererV2)
22
23extern AliRun *gAlice;
24
25AliITSclustererV2::AliITSclustererV2(const AliITSgeom *geom) {
26 //------------------------------------------------------------
27 // Standard constructor
28 //------------------------------------------------------------
29 AliITSgeom *g=(AliITSgeom*)geom;
30
31 fEvent=0;
32 fI=0;
33
34 Int_t mmax=geom->GetIndexMax();
35 if (mmax>2200) {cerr<<"Too many ITS subdetectors !\n"; exit(1);}
36 Int_t m;
37 for (m=0; m<mmax; m++) {
38 Int_t lay,lad,det; g->GetModuleId(m,lay,lad,det);
39 Float_t x,y,z; g->GetTrans(lay,lad,det,x,y,z);
40 Double_t rot[9]; g->GetRotMatrix(lay,lad,det,rot);
41 fYshift[m] = x*rot[0] + y*rot[1];
42 fZshift[m] = (Double_t)z;
43 fNdet[m] = (lad-1)*g->GetNdetectors(lay) + (det-1);
44 }
45
46 //SPD geometry
47 fLastSPD1=g->GetModuleIndex(2,1,1)-1;
48 fNySPD=256; fNzSPD=160;
49 fYpitchSPD=0.0050;
50 fZ1pitchSPD=0.0425; fZ2pitchSPD=0.0625;
51 fHwSPD=0.64; fHlSPD=3.48;
52 fYSPD[0]=0.5*fYpitchSPD;
53 for (m=1; m<fNySPD; m++) fYSPD[m]=fYSPD[m-1]+fYpitchSPD;
54 fZSPD[0]=fZ1pitchSPD;
55 for (m=1; m<fNzSPD; m++) {
56 Double_t dz=fZ1pitchSPD;
57 if (m==31 || m==32 || m==63 || m==64 || m==95 || m==96 ||
58 m==127 || m==128) dz=fZ2pitchSPD;
59 fZSPD[m]=fZSPD[m-1]+dz;
60 }
61 for (m=0; m<fNzSPD; m++) {
62 Double_t dz=0.5*fZ1pitchSPD;
63 if (m==31 || m==32 || m==63 || m==64 || m==95 || m==96 ||
64 m==127 || m==128) dz=0.5*fZ2pitchSPD;
65 fZSPD[m]-=dz;
66 }
67
68 //SDD geometry
69 fNySDD=256; fNzSDD=256;
70 fYpitchSDD=0.01825;
71 fZpitchSDD=0.02940;
72 fHwSDD=3.5085; fHlSDD=3.7632;
73 fYoffSDD=0.0425;
74
75 //SSD geometry
76 fLastSSD1=g->GetModuleIndex(6,1,1)-1;
77 fYpitchSSD=0.0095;
78 fHwSSD=3.65;
79 fHlSSD=2.00;
80 fTanP=0.0275;
81 fTanN=0.0075;
82}
83
84void AliITSclustererV2::Digits2Clusters(const TFile *in, TFile *out) {
85 //------------------------------------------------------------
86 // This function creates ITS clusters
87 //------------------------------------------------------------
88 Int_t ncl=0;
89 TDirectory *savedir=gDirectory;
90
91 if (!out->IsOpen()) {
92 cerr<<"AliITSclustererV2::Digits2Clusters(): output file not open !\n";
93 return;
94 }
95
96 Char_t name[100];
97 sprintf(name,"TreeD%d",fEvent);
98
99 //TTree *dTree=(TTree*)((TFile*)in)->Get(name);
100 TTree *dTree=gAlice->TreeD();
101 if (!dTree) {
102 cerr<<"Input tree "<<name<<" not found !\n";
103 return;
104 }
105
106 TClonesArray *digitsSPD=new TClonesArray("AliITSdigitSPD",3000);
107 dTree->SetBranchAddress("ITSDigitsSPD",&digitsSPD);
108 TClonesArray *digitsSDD=new TClonesArray("AliITSdigitSDD",3000);
109 dTree->SetBranchAddress("ITSDigitsSDD",&digitsSDD);
110 TClonesArray *digitsSSD=new TClonesArray("AliITSdigitSSD",3000);
111 dTree->SetBranchAddress("ITSDigitsSSD",&digitsSSD);
112
113 Int_t mmax=(Int_t)dTree->GetEntries();
114
115 out->cd();
116
117 sprintf(name,"TreeC_ITS_%d",fEvent);
118 TTree cTree(name,"ITS clusters V2");
119 TClonesArray *clusters=new TClonesArray("AliITSclusterV2",1000);
120 cTree.Branch("Clusters",&clusters);
121
122 for (fI=0; fI<mmax; fI++) {
123 dTree->GetEvent(fI);
124
125 if (digitsSPD->GetEntriesFast()!=0)
126 FindClustersSPD(digitsSPD,clusters);
127 else if(digitsSDD->GetEntriesFast()!=0)
128 FindClustersSDD(digitsSDD,clusters);
129 else if(digitsSSD->GetEntriesFast()!=0)
130 FindClustersSSD(digitsSSD,clusters);
131
132 ncl+=clusters->GetEntriesFast();
133
134 cTree.Fill();
135
136 digitsSPD->Clear();
137 digitsSDD->Clear();
138 digitsSSD->Clear();
139 clusters->Clear();
140 }
141 cTree.Write();
142
143 delete clusters;
144
145 delete digitsSPD;
146 delete digitsSDD;
147 delete digitsSSD;
148
149 //delete dTree;
150
151 cerr<<"Number of found clusters : "<<ncl<<endl;
152
153 savedir->cd();
154}
155
156//**** Fast clusters *******************************
157#include "TParticle.h"
158
159#include "AliITS.h"
160#include "AliITSmodule.h"
161#include "AliITSRecPoint.h"
162#include "AliITSsimulationFastPoints.h"
163#include "AliITSRecPoint.h"
164
165static void CheckLabels(Int_t lab[3]) {
166 //------------------------------------------------------------
167 // Tries to find mother's labels
168 //------------------------------------------------------------
169 Int_t label=lab[0];
170 if (label>=0) {
171 TParticle *part=(TParticle*)gAlice->Particle(label);
172 label=-3;
173 while (part->P() < 0.005) {
174 Int_t m=part->GetFirstMother();
175 if (m<0) {cerr<<"Primary momentum: "<<part->P()<<endl; break;}
1cca57bf 176 if (part->GetStatusCode()>0) {
177 cerr<<"Primary momentum: "<<part->P()<<endl; break;
178 }
7f4044f0 179 label=m;
180 part=(TParticle*)gAlice->Particle(label);
181 }
182 if (lab[1]<0) lab[1]=label;
183 else if (lab[2]<0) lab[2]=label;
184 else ;//cerr<<"CheckLabels : No empty labels !\n";
185 }
186}
187
1cca57bf 188void AliITSclustererV2::RecPoints2Clusters
189(const TClonesArray *points, Int_t idx, TClonesArray *clusters) {
7f4044f0 190 //------------------------------------------------------------
1cca57bf 191 // Conversion AliITSRecPoint -> AliITSclusterV2 for the ITS
192 // subdetector indexed by idx
7f4044f0 193 //------------------------------------------------------------
194 TClonesArray &cl=*clusters;
195 Int_t ncl=points->GetEntriesFast();
196 for (Int_t i=0; i<ncl; i++) {
197 AliITSRecPoint *p = (AliITSRecPoint *)points->UncheckedAt(i);
198 Float_t lp[5];
1cca57bf 199 lp[0]=-p->GetX()-fYshift[idx]; if (idx<=fLastSPD1) lp[0]*=-1; //SPD1
200 lp[1]=p->GetZ()+fZshift[idx];
201 lp[2]=p->GetSigmaX2();
202 lp[3]=p->GetSigmaZ2();
203 lp[4]=p->GetQ();
7f4044f0 204 Int_t lab[4];
205 lab[0]=p->GetLabel(0); lab[1]=p->GetLabel(1); lab[2]=p->GetLabel(2);
1cca57bf 206 lab[3]=fNdet[idx];
7f4044f0 207 CheckLabels(lab);
208 new (cl[i]) AliITSclusterV2(lab,lp);
209 }
210}
211
212void AliITSclustererV2::Hits2Clusters(const TFile *in, TFile *out) {
213 //------------------------------------------------------------
214 // This function creates ITS clusters
215 //------------------------------------------------------------
216 TDirectory *savedir=gDirectory;
217
218 if (!out->IsOpen()) {
219 cerr<<"AliITSclustererV2::Hits2Clusters: output file not open !\n";
220 return;
221 }
222
223 if (!gAlice) {
224 cerr<<"AliITSclustererV2::Hits2Clusters : gAlice==0 !\n";
225 return;
226 }
227
228 AliITS *its = (AliITS*)gAlice->GetModule("ITS");
229 if (!its) {
230 cerr<<"AliITSclustererV2::Hits2Clusters : Can't find the ITS !\n";
231 return;
232 }
233 AliITSgeom *geom=its->GetITSgeom();
234 Int_t mmax=geom->GetIndexMax();
235
236 its->InitModules(-1,mmax);
237 its->FillModules(gAlice->TreeH(),0);
238
239 out->cd();
240
241 Char_t name[100];
242 sprintf(name,"TreeC_ITS_%d",fEvent);
243 TTree cTree(name,"ITS clusters V2");
244 TClonesArray *clusters=new TClonesArray("AliITSclusterV2",1000);
245 cTree.Branch("Clusters",&clusters);
246
247 static TClonesArray *points=its->RecPoints();
248 AliITSsimulationFastPoints sim;
249 Int_t ncl=0;
250 for (Int_t m=0; m<mmax; m++) {
251 AliITSmodule *mod=its->GetModule(m);
252 sim.CreateFastRecPoints(mod,m,gRandom);
253
1cca57bf 254 RecPoints2Clusters(points, m, clusters);
7f4044f0 255 its->ResetRecPoints();
256
257 ncl+=clusters->GetEntriesFast();
258 cTree.Fill();
259 clusters->Clear();
260 }
261 cTree.Write();
262
263 cerr<<"Number of found fast clusters : "<<ncl<<endl;
264
265 delete clusters;
266
267 savedir->cd();
268}
269
270//***********************************
271
272#ifndef V1
273
274void AliITSclustererV2::
275FindCluster(Int_t k,Int_t maxz,AliBin *bins,Int_t &n,Int_t *idx) {
276 //------------------------------------------------------------
277 // returns an array of indices of digits belonging to the cluster
278 // (needed when the segmentation is not regular)
279 //------------------------------------------------------------
280 if (n<200) idx[n++]=bins[k].GetIndex();
281 bins[k].Use();
282
283 if (bins[k-maxz].IsNotUsed()) FindCluster(k-maxz,maxz,bins,n,idx);
284 if (bins[k-1 ].IsNotUsed()) FindCluster(k-1 ,maxz,bins,n,idx);
285 if (bins[k+maxz].IsNotUsed()) FindCluster(k+maxz,maxz,bins,n,idx);
286 if (bins[k+1 ].IsNotUsed()) FindCluster(k+1 ,maxz,bins,n,idx);
287 /*
288 if (bins[k-maxz-1].IsNotUsed()) FindCluster(k-maxz-1,maxz,bins,n,idx);
289 if (bins[k-maxz+1].IsNotUsed()) FindCluster(k-maxz+1,maxz,bins,n,idx);
290 if (bins[k+maxz-1].IsNotUsed()) FindCluster(k+maxz-1,maxz,bins,n,idx);
291 if (bins[k+maxz+1].IsNotUsed()) FindCluster(k+maxz+1,maxz,bins,n,idx);
292 */
293}
294
295void AliITSclustererV2::
296FindClustersSPD(const TClonesArray *digits, TClonesArray *clusters) {
297 //------------------------------------------------------------
298 // Actual SPD cluster finder
299 //------------------------------------------------------------
300 const Int_t kMAXBIN=(fNzSPD+2)*(fNySPD+2);
301
302 Int_t ndigits=digits->GetEntriesFast();
303 AliBin *bins=new AliBin[kMAXBIN];
304
305 Int_t k;
306 AliITSdigitSPD *d=0;
307 for (k=0; k<ndigits; k++) {
308 d=(AliITSdigitSPD*)digits->UncheckedAt(k);
309 Int_t i=d->GetCoord2()+1; //y
310 Int_t j=d->GetCoord1()+1;
311 bins[i*fNzSPD+j].SetIndex(k);
312 bins[i*fNzSPD+j].SetMask(1);
313 }
314
315 Int_t n=0; TClonesArray &cl=*clusters;
316 for (k=0; k<kMAXBIN; k++) {
317 if (!bins[k].IsNotUsed()) continue;
318 Int_t ni=0, idx[200];
319 FindCluster(k,fNzSPD,bins,ni,idx);
320 if (ni==200) {cerr<<"SPD: Too big cluster !\n"; continue;}
321
322 Int_t lab[4];
323 lab[0]=-2;
324 lab[1]=-2;
325 lab[2]=-2;
326 lab[3]=fNdet[fI];
327
328 d=(AliITSdigitSPD*)digits->UncheckedAt(idx[0]);
329 Int_t ymin=d->GetCoord2(),ymax=ymin;
330 Int_t zmin=d->GetCoord1(),zmax=zmin;
331 Float_t y=0.,z=0.,q=0.;
332 for (Int_t l=0; l<ni; l++) {
333 d=(AliITSdigitSPD*)digits->UncheckedAt(idx[l]);
334
335 if (ymin > d->GetCoord2()) ymin=d->GetCoord2();
336 if (ymax < d->GetCoord2()) ymax=d->GetCoord2();
337 if (zmin > d->GetCoord1()) zmin=d->GetCoord1();
338 if (zmax < d->GetCoord1()) zmax=d->GetCoord1();
339
340 Int_t lab0=(d->GetTracks())[0];
341 if (lab0>=0) {
342 if (lab[0]<0) {
343 lab[0]=lab0;
344 } else if (lab[1]<0) {
345 if (lab0!=lab[0]) lab[1]=lab0;
346 } else if (lab[2]<0) {
347 if (lab0!=lab[0])
348 if (lab0!=lab[1]) lab[2]=lab0;
349 }
350 }
351 Float_t qq=d->GetSignal();
352 y+=qq*fYSPD[d->GetCoord2()]; z+=qq*fZSPD[d->GetCoord1()]; q+=qq;
353 }
354 y/=q; z/=q;
355 y-=fHwSPD; z-=fHlSPD;
356
357 Float_t lp[5];
358 lp[0]=-y-fYshift[fI]; if (fI<=fLastSPD1) lp[0]=-lp[0];
359 lp[1]= z+fZshift[fI];
360 lp[2]= fYpitchSPD*fYpitchSPD/12.;
361 lp[3]= fZ1pitchSPD*fZ1pitchSPD/12.;
362 //lp[4]= q;
363 lp[4]= (zmax-zmin+1)*100 + (ymax-ymin+1);
364
365 //CheckLabels(lab);
366 new (cl[n]) AliITSclusterV2(lab,lp); n++;
367 }
368
369 delete bins;
370}
371
372Bool_t AliITSclustererV2::IsMaximum(Int_t k,Int_t max,const AliBin *bins) {
373 //------------------------------------------------------------
374 //is this a local maximum ?
375 //------------------------------------------------------------
376 UShort_t q=bins[k].GetQ();
377 if (q==1023) return kFALSE;
378 if (bins[k-max].GetQ() > q) return kFALSE;
379 if (bins[k-1 ].GetQ() > q) return kFALSE;
380 if (bins[k+max].GetQ() > q) return kFALSE;
381 if (bins[k+1 ].GetQ() > q) return kFALSE;
382 if (bins[k-max-1].GetQ() > q) return kFALSE;
383 if (bins[k+max-1].GetQ() > q) return kFALSE;
384 if (bins[k+max+1].GetQ() > q) return kFALSE;
385 if (bins[k-max+1].GetQ() > q) return kFALSE;
386 return kTRUE;
387}
388
389void AliITSclustererV2::
390FindPeaks(Int_t k,Int_t max,AliBin *b,Int_t *idx,UInt_t *msk,Int_t& n) {
391 //------------------------------------------------------------
392 //find local maxima
393 //------------------------------------------------------------
394 if (n<31)
395 if (IsMaximum(k,max,b)) {
396 idx[n]=k; msk[n]=(2<<n);
397 n++;
398 }
399 b[k].SetMask(0);
400 if (b[k-max].GetMask()&1) FindPeaks(k-max,max,b,idx,msk,n);
401 if (b[k-1 ].GetMask()&1) FindPeaks(k-1 ,max,b,idx,msk,n);
402 if (b[k+max].GetMask()&1) FindPeaks(k+max,max,b,idx,msk,n);
403 if (b[k+1 ].GetMask()&1) FindPeaks(k+1 ,max,b,idx,msk,n);
404}
405
406void AliITSclustererV2::
407MarkPeak(Int_t k, Int_t max, AliBin *bins, UInt_t m) {
408 //------------------------------------------------------------
409 //mark this peak
410 //------------------------------------------------------------
411 UShort_t q=bins[k].GetQ();
412
413 bins[k].SetMask(bins[k].GetMask()|m);
414
415 if (bins[k-max].GetQ() <= q)
416 if ((bins[k-max].GetMask()&m) == 0) MarkPeak(k-max,max,bins,m);
417 if (bins[k-1 ].GetQ() <= q)
418 if ((bins[k-1 ].GetMask()&m) == 0) MarkPeak(k-1 ,max,bins,m);
419 if (bins[k+max].GetQ() <= q)
420 if ((bins[k+max].GetMask()&m) == 0) MarkPeak(k+max,max,bins,m);
421 if (bins[k+1 ].GetQ() <= q)
422 if ((bins[k+1 ].GetMask()&m) == 0) MarkPeak(k+1 ,max,bins,m);
423}
424
425void AliITSclustererV2::
426MakeCluster(Int_t k,Int_t max,AliBin *bins,UInt_t m,AliITSclusterV2 &c) {
427 //------------------------------------------------------------
428 //make cluster using digits of this peak
429 //------------------------------------------------------------
430 Float_t q=(Float_t)bins[k].GetQ();
431 Int_t i=k/max, j=k-i*max;
432
433 c.SetQ(c.GetQ()+q);
434 c.SetY(c.GetY()+i*q);
435 c.SetZ(c.GetZ()+j*q);
436 c.SetSigmaY2(c.GetSigmaY2()+i*i*q);
437 c.SetSigmaZ2(c.GetSigmaZ2()+j*j*q);
438
439 bins[k].SetMask(0xFFFFFFFE);
440
441 if (bins[k-max].GetMask() == m) MakeCluster(k-max,max,bins,m,c);
442 if (bins[k-1 ].GetMask() == m) MakeCluster(k-1 ,max,bins,m,c);
443 if (bins[k+max].GetMask() == m) MakeCluster(k+max,max,bins,m,c);
444 if (bins[k+1 ].GetMask() == m) MakeCluster(k+1 ,max,bins,m,c);
445}
446
447void AliITSclustererV2::
448FindClustersSDD(const TClonesArray *digits, TClonesArray *clusters) {
449 //------------------------------------------------------------
450 // Actual SDD cluster finder
451 //------------------------------------------------------------
452 const Int_t kMAXBIN=(fNzSDD+2)*(fNySDD+2);
453
454 AliBin *bins[2];
455 bins[0]=new AliBin[kMAXBIN];
456 bins[1]=new AliBin[kMAXBIN];
457
458 AliITSdigitSDD *d=0;
459 Int_t i, ndigits=digits->GetEntriesFast();
460 for (i=0; i<ndigits; i++) {
461 d=(AliITSdigitSDD*)digits->UncheckedAt(i);
462 Int_t y=d->GetCoord2()+1; //y
463 Int_t z=d->GetCoord1()+1; //z
464 Int_t q=d->GetSignal();
465 if (z <= fNzSDD) {
466 bins[0][y*fNzSDD+z].SetQ(q);
467 bins[0][y*fNzSDD+z].SetMask(1);
468 bins[0][y*fNzSDD+z].SetIndex(i);
469 } else {
470 z-=fNzSDD;
471 bins[1][y*fNzSDD+z].SetQ(q);
472 bins[1][y*fNzSDD+z].SetMask(1);
473 bins[1][y*fNzSDD+z].SetIndex(i);
474 }
475 }
476
477 Int_t ncl=0; TClonesArray &cl=*clusters;
478 for (Int_t s=0; s<2; s++)
479 for (i=0; i<kMAXBIN; i++) {
480 if (bins[s][i].IsUsed()) continue;
481 Int_t idx[32]; UInt_t msk[32]; Int_t npeaks=0;
482 FindPeaks(i, fNzSDD, bins[s], idx, msk, npeaks);
483
484 if (npeaks>30) continue;
485
486 Int_t k,l;
487 for (k=0; k<npeaks-1; k++){//mark adjacent peaks
488 if (idx[k] < 0) continue; //this peak is already removed
489 for (l=k+1; l<npeaks; l++) {
490 if (idx[l] < 0) continue; //this peak is already removed
491 Int_t ki=idx[k]/fNzSDD, kj=idx[k] - ki*fNzSDD;
492 Int_t li=idx[l]/fNzSDD, lj=idx[l] - li*fNzSDD;
493 Int_t di=TMath::Abs(ki - li);
494 Int_t dj=TMath::Abs(kj - lj);
495 if (di>1 || dj>1) continue;
496 if (bins[s][idx[k]].GetQ() > bins[s][idx[l]].GetQ()) {
497 msk[l]=msk[k];
498 idx[l]*=-1;
499 } else {
500 msk[k]=msk[l];
501 idx[k]*=-1;
502 break;
503 }
504 }
505 }
506
507 for (k=0; k<npeaks; k++) {
508 MarkPeak(TMath::Abs(idx[k]), fNzSDD, bins[s], msk[k]);
509 }
510
511 for (k=0; k<npeaks; k++) {
512 if (idx[k] < 0) continue; //removed peak
513 AliITSclusterV2 c;
514 MakeCluster(idx[k], fNzSDD, bins[s], msk[k], c);
515
516 //if (c.GetQ() < 200) continue; //noise cluster
517
518 /*
519 Float_t s2 = c.GetSigmaY2()/c.GetQ() - c.GetY()*c.GetY();
520 Float_t w=par->GetPadPitchWidth(sec);
521 c.SetSigmaY2((s2 + 1./12.)*w*w);
522 if (s2 != 0.) {
523 c.SetSigmaY2(c.GetSigmaY2()*0.108);
524 if (sec<par->GetNInnerSector()) c.SetSigmaY2(c.GetSigmaY2()*2.07);
525 }
526
527 s2 = c.GetSigmaZ2()/c.GetQ() - c.GetZ()*c.GetZ();
528 w=par->GetZWidth();
529 c.SetSigmaZ2((s2 + 1./12.)*w*w);
530 if (s2 != 0.) {
531 c.SetSigmaZ2(c.GetSigmaZ2()*0.169);
532 if (sec<par->GetNInnerSector()) c.SetSigmaZ2(c.GetSigmaZ2()*1.77);
533 }
534 */
535
536 c.SetSigmaY2(0.0030*0.0030);
537 c.SetSigmaZ2(0.0020*0.0020);
538 c.SetDetectorIndex(fNdet[fI]);
539
540 Float_t y=c.GetY(),z=c.GetZ(), q=c.GetQ();
541 y/=q; z/=q;
542
543 y=(y-0.5)*fYpitchSDD;
544 y-=fHwSDD;
545 y-=fYoffSDD; //delay ?
546 if (s) y=-y;
547
548 z=(z-0.5)*fZpitchSDD;
549 z-=fHlSDD;
550
551 y=-y-fYshift[fI];
552 z= z+fZshift[fI];
553 c.SetY(y);
554 c.SetZ(z);
555
556 c.SetQ(q/20.); //to be consistent with the SSD charges
557
558 AliBin *b=&bins[s][idx[k]];
559 d=(AliITSdigitSDD*)digits->UncheckedAt(b->GetIndex());
560 Int_t l0=(d->GetTracks())[0];
561 if (l0<0) {
562 b=&bins[s][idx[k]-1];
563 if (b->GetQ()>0) {
564 d=(AliITSdigitSDD*)digits->UncheckedAt(b->GetIndex());
565 l0=(d->GetTracks())[0];
566 }
567 }
568 if (l0<0) {
569 b=&bins[s][idx[k]+1];
570 if (b->GetQ()>0) {
571 d=(AliITSdigitSDD*)digits->UncheckedAt(b->GetIndex());
572 l0=(d->GetTracks())[0];
573 }
574 }
575 if (l0<0) {
576 b=&bins[s][idx[k]-fNzSDD];
577 if (b->GetQ()>0) {
578 d=(AliITSdigitSDD*)digits->UncheckedAt(b->GetIndex());
579 l0=(d->GetTracks())[0];
580 }
581 }
582 if (l0<0) {
583 b=&bins[s][idx[k]+fNzSDD];
584 if (b->GetQ()>0) {
585 d=(AliITSdigitSDD*)digits->UncheckedAt(b->GetIndex());
586 l0=(d->GetTracks())[0];
587 }
588 }
589
590 if (l0<0) {
591 b=&bins[s][idx[k]+fNzSDD+1];
592 if (b->GetQ()>0) {
593 d=(AliITSdigitSDD*)digits->UncheckedAt(b->GetIndex());
594 l0=(d->GetTracks())[0];
595 }
596 }
597 if (l0<0) {
598 b=&bins[s][idx[k]+fNzSDD-1];
599 if (b->GetQ()>0) {
600 d=(AliITSdigitSDD*)digits->UncheckedAt(b->GetIndex());
601 l0=(d->GetTracks())[0];
602 }
603 }
604 if (l0<0) {
605 b=&bins[s][idx[k]-fNzSDD+1];
606 if (b->GetQ()>0) {
607 d=(AliITSdigitSDD*)digits->UncheckedAt(b->GetIndex());
608 l0=(d->GetTracks())[0];
609 }
610 }
611 if (l0<0) {
612 b=&bins[s][idx[k]-fNzSDD-1];
613 if (b->GetQ()>0) {
614 d=(AliITSdigitSDD*)digits->UncheckedAt(b->GetIndex());
615 l0=(d->GetTracks())[0];
616 }
617 }
618
619 {
620 Int_t lab[3];
621 lab[0]=(d->GetTracks())[0];
622 lab[1]=(d->GetTracks())[1];
623 lab[2]=(d->GetTracks())[2];
624 //CheckLabels(lab);
625 c.SetLabel(lab[0],0);
626 c.SetLabel(lab[1],1);
627 c.SetLabel(lab[2],2);
628 }
629
630 new (cl[ncl]) AliITSclusterV2(c); ncl++;
631 }
632 }
633
634 delete[] bins[0];
635 delete[] bins[1];
636}
637
638void AliITSclustererV2::
639FindClustersSSD(const TClonesArray *digits, TClonesArray *clusters) {
640 //------------------------------------------------------------
641 // Actual SSD cluster finder
642 //------------------------------------------------------------
643 Int_t smax=digits->GetEntriesFast();
644 if (smax==0) return;
645
646 const Int_t MAX=1000;
647 Int_t np=0, nn=0;
648 Ali1Dcluster pos[MAX], neg[MAX];
649 Float_t y=0., q=0., qmax=0.;
650 Int_t lab[4]={-2,-2,-2,-2};
651
652 TClonesArray &cl=*clusters;
653
654 AliITSdigitSSD *d=(AliITSdigitSSD*)digits->UncheckedAt(0);
655 q += d->GetSignal();
656 y += d->GetCoord2()*d->GetSignal();
657 qmax=d->GetSignal();
658 lab[0]=d->GetTrack(0); lab[1]=d->GetTrack(1); lab[2]=d->GetTrack(2);
659 Int_t curr=d->GetCoord2();
660 Int_t flag=d->GetCoord1();
661 Int_t *n=&nn;
662 Ali1Dcluster *c=neg;
663 Int_t nd=0;
664 for (Int_t s=1; s<smax; s++) {
665 d=(AliITSdigitSSD*)digits->UncheckedAt(s);
666 Int_t strip=d->GetCoord2();
667 if ((strip-curr) > 1 || flag!=d->GetCoord1()) {
668 c[*n].SetY(y/q);
669 c[*n].SetQ(q);
670 c[*n].SetNd(nd);
671 c[*n].SetLabels(lab);
672 //Split suspiciously big cluster
673 if (nd>3) {
674 c[*n].SetY(y/q-0.5*nd);
675 c[*n].SetQ(0.5*q);
676 (*n)++;
677 if (*n==MAX) {
678 cerr<<
679 "AliITSclustererV2::FindClustersSSD: Too many 1D clusters !\n";
680 return;
681 }
682 c[*n].SetY(y/q+0.5*nd);
683 c[*n].SetQ(0.5*q);
684 c[*n].SetNd(nd);
685 c[*n].SetLabels(lab);
686 }
687 (*n)++;
688 if (*n==MAX) {
689 cerr<<"AliITSclustererV2::FindClustersSSD: Too many 1D clusters !\n";
690 return;
691 }
692 y=q=qmax=0.;
693 nd=0;
694 lab[0]=lab[1]=lab[2]=-2;
695 if (flag!=d->GetCoord1()) { n=&np; c=pos; }
696 }
697 flag=d->GetCoord1();
698 q += d->GetSignal();
699 y += d->GetCoord2()*d->GetSignal();
700 nd++;
701 if (d->GetSignal()>qmax) {
702 qmax=d->GetSignal();
703 lab[0]=d->GetTrack(0); lab[1]=d->GetTrack(1); lab[2]=d->GetTrack(2);
704 }
705 curr=strip;
706 }
707 c[*n].SetY(y/q);
708 c[*n].SetQ(q);
709 c[*n].SetNd(nd);
710 c[*n].SetLabels(lab);
711 //Split suspiciously big cluster
712 if (nd>3) {
713 c[*n].SetY(y/q-0.5*nd);
714 c[*n].SetQ(0.5*q);
715 (*n)++;
716 if (*n==MAX) {
717 cerr<<"AliITSclustererV2::FindClustersSSD: Too many 1D clusters !\n";
718 return;
719 }
720 c[*n].SetY(y/q+0.5*nd);
721 c[*n].SetQ(0.5*q);
722 c[*n].SetNd(nd);
723 c[*n].SetLabels(lab);
724 }
725 (*n)++;
726 if (*n==MAX) {
727 cerr<<"AliITSclustererV2::FindClustersSSD: Too many 1D clusters !\n";
728 return;
729 }
730
731 Float_t tanp=fTanP, tann=fTanN;
732 if (fI>fLastSSD1) {tann=fTanP; tanp=fTanN;}
733
734 Int_t idet=fNdet[fI];
735 Int_t ncl=0;
736 for (Int_t i=0; i<np; i++) {
737 //Float_t dq_min=1.e+33;
738 Float_t ybest=1000,zbest=1000,qbest=0;
739 Float_t yp=pos[i].GetY()*fYpitchSSD;
740 for (Int_t j=0; j<nn; j++) {
741 //if (pos[i].fTracks[0] != neg[j].fTracks[0]) continue;
742 Float_t yn=neg[j].GetY()*fYpitchSSD;
743 Float_t zt=(2*fHlSSD*tanp + yp - yn)/(tann+tanp);
744 Float_t yt=yn + tann*zt;
745 zt-=fHlSSD; yt-=fHwSSD;
746 if (TMath::Abs(yt)<fHwSSD+0.01)
747 if (TMath::Abs(zt)<fHlSSD+0.01) {
748 //if (TMath::Abs(pos[i].GetQ()-neg[j].GetQ())<dq_min) {
749 //dq_min=TMath::Abs(pos[i].GetQ()-neg[j].GetQ());
750 ybest=yt; zbest=zt;
751 qbest=0.5*(pos[i].GetQ()+neg[j].GetQ());
752
753 lab[0]=pos[i].GetLabel(0);
754 lab[1]=pos[i].GetLabel(1);
755 lab[2]=neg[i].GetLabel(0);
756 lab[3]=(((i<<10) + j)<<10) + idet; // pos|neg|det
757 Float_t lp[5];
758 lp[0]=-ybest-fYshift[fI];
759 lp[1]= zbest+fZshift[fI];
760 lp[2]=0.0025*0.0025; //SigmaY2
761 lp[3]=0.110*0.110; //SigmaZ2
762 if (pos[i].GetNd()+neg[j].GetNd() > 4) {
763 lp[2]*=9;
764 lp[3]*=9;
765 }
766 lp[4]=qbest; //Q
767
768 //CheckLabels(lab);
769 new (cl[ncl]) AliITSclusterV2(lab,lp); ncl++;
770 }
771 }
772 /*
773 if (ybest<100) {
774 lab[3]=idet;
775 Float_t lp[5];
776 lp[0]=-ybest-fYshift[fI];
777 lp[1]= zbest+fZshift[fI];
778 lp[2]=0.002*0.002; //SigmaY2
779 lp[3]=0.080*0.080; //SigmaZ2
780 lp[4]=qbest; //Q
781 //
782 new (cl[ncl]) AliITSclusterV2(lab,lp); ncl++;
783 }
784 */
785 }
786}
787
788#else //V1
789
790#include "AliITSDetType.h"
791
792#include "AliITSsegmentationSPD.h"
793#include "AliITSClusterFinderSPD.h"
794
795#include "AliITSresponseSDD.h"
796#include "AliITSsegmentationSDD.h"
797#include "AliITSClusterFinderSDD.h"
798
799#include "AliITSsegmentationSSD.h"
800#include "AliITSClusterFinderSSD.h"
801
802
803void AliITSclustererV2::
804FindClustersSPD(const TClonesArray *digits, TClonesArray *clusters) {
805 //------------------------------------------------------------
806 // Actual SPD cluster finding based on AliITSClusterFinderSPD
807 //------------------------------------------------------------
808 static AliITS *its=(AliITS*)gAlice->GetModule("ITS");
809 static TClonesArray *points=its->RecPoints();
810 static AliITSsegmentationSPD *seg=
811 (AliITSsegmentationSPD *)its->DetType(0)->GetSegmentationModel();
812 static AliITSClusterFinderSPD cf(seg, (TClonesArray*)digits, points);
813
814 cf.FindRawClusters(fI);
1cca57bf 815 RecPoints2Clusters(points, fI, clusters);
7f4044f0 816 its->ResetRecPoints();
817
818}
819
820void AliITSclustererV2::
821FindClustersSDD(const TClonesArray *digits, TClonesArray *clusters) {
822 //------------------------------------------------------------
823 // Actual SDD cluster finding based on AliITSClusterFinderSDD
824 //------------------------------------------------------------
825 static AliITS *its=(AliITS*)gAlice->GetModule("ITS");
826 static TClonesArray *points=its->RecPoints();
827 static AliITSresponseSDD *resp=
828 (AliITSresponseSDD *)its->DetType(1)->GetResponseModel();
829 static AliITSsegmentationSDD *seg=
830 (AliITSsegmentationSDD *)its->DetType(1)->GetSegmentationModel();
831 static AliITSClusterFinderSDD
832 cf(seg,resp,(TClonesArray*)digits,its->ClustersAddress(1));
833
834 cf.FindRawClusters(fI);
1cca57bf 835 Int_t nc=points->GetEntriesFast();
836 while (nc--) { //To be consistent with the SSD cluster charges
837 AliITSRecPoint *p=(AliITSRecPoint*)points->UncheckedAt(nc);
838 p->SetQ(p->GetQ/12.);
839 }
840 RecPoints2Clusters(points, fI, clusters);
7f4044f0 841 its->ResetClusters(1);
842 its->ResetRecPoints();
843
844}
845
846void AliITSclustererV2::
847FindClustersSSD(const TClonesArray *digits, TClonesArray *clusters) {
848 //------------------------------------------------------------
849 // Actual SSD cluster finding based on AliITSClusterFinderSSD
850 //------------------------------------------------------------
851 static AliITS *its=(AliITS*)gAlice->GetModule("ITS");
852 static TClonesArray *points=its->RecPoints();
853 static AliITSsegmentationSSD *seg=
854 (AliITSsegmentationSSD *)its->DetType(2)->GetSegmentationModel();
855 static AliITSClusterFinderSSD cf(seg,(TClonesArray*)digits);
856
857 cf.FindRawClusters(fI);
1cca57bf 858 RecPoints2Clusters(points, fI, clusters);
7f4044f0 859 its->ResetRecPoints();
860
861}
862
863#endif