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