]>
Commit | Line | Data |
---|---|---|
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 | ||
21 | ClassImp(AliITSclustererV2) | |
22 | ||
23 | extern AliRun *gAlice; | |
24 | ||
25 | AliITSclustererV2::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 | ||
84 | void 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 | ||
165 | static 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;} | |
176 | label=m; | |
177 | part=(TParticle*)gAlice->Particle(label); | |
178 | } | |
179 | if (lab[1]<0) lab[1]=label; | |
180 | else if (lab[2]<0) lab[2]=label; | |
181 | else ;//cerr<<"CheckLabels : No empty labels !\n"; | |
182 | } | |
183 | } | |
184 | ||
185 | static void RecPoints2Clusters(const TClonesArray *points, | |
186 | Double_t ys, Double_t zs, Int_t idx, Int_t sgn, TClonesArray *clusters) { | |
187 | //------------------------------------------------------------ | |
188 | // Conversion AliITSRecPoint -> AliITSclusterV2 | |
189 | //------------------------------------------------------------ | |
190 | TClonesArray &cl=*clusters; | |
191 | Int_t ncl=points->GetEntriesFast(); | |
192 | for (Int_t i=0; i<ncl; i++) { | |
193 | AliITSRecPoint *p = (AliITSRecPoint *)points->UncheckedAt(i); | |
194 | Float_t lp[5]; | |
195 | lp[0]=-p->GetX()-ys; lp[0]*=sgn; //SPD1 | |
196 | lp[1]=p->GetZ()+zs; | |
197 | lp[2]=p->GetSigmaX2()*0.5; | |
198 | lp[3]=p->GetSigmaZ2()*0.5; | |
199 | lp[4]=p->GetQ()/12.; //to be roughly consistent with SSD | |
200 | Int_t lab[4]; | |
201 | lab[0]=p->GetLabel(0); lab[1]=p->GetLabel(1); lab[2]=p->GetLabel(2); | |
202 | lab[3]=idx; | |
203 | CheckLabels(lab); | |
204 | new (cl[i]) AliITSclusterV2(lab,lp); | |
205 | } | |
206 | } | |
207 | ||
208 | void AliITSclustererV2::Hits2Clusters(const TFile *in, TFile *out) { | |
209 | //------------------------------------------------------------ | |
210 | // This function creates ITS clusters | |
211 | //------------------------------------------------------------ | |
212 | TDirectory *savedir=gDirectory; | |
213 | ||
214 | if (!out->IsOpen()) { | |
215 | cerr<<"AliITSclustererV2::Hits2Clusters: output file not open !\n"; | |
216 | return; | |
217 | } | |
218 | ||
219 | if (!gAlice) { | |
220 | cerr<<"AliITSclustererV2::Hits2Clusters : gAlice==0 !\n"; | |
221 | return; | |
222 | } | |
223 | ||
224 | AliITS *its = (AliITS*)gAlice->GetModule("ITS"); | |
225 | if (!its) { | |
226 | cerr<<"AliITSclustererV2::Hits2Clusters : Can't find the ITS !\n"; | |
227 | return; | |
228 | } | |
229 | AliITSgeom *geom=its->GetITSgeom(); | |
230 | Int_t mmax=geom->GetIndexMax(); | |
231 | ||
232 | its->InitModules(-1,mmax); | |
233 | its->FillModules(gAlice->TreeH(),0); | |
234 | ||
235 | out->cd(); | |
236 | ||
237 | Char_t name[100]; | |
238 | sprintf(name,"TreeC_ITS_%d",fEvent); | |
239 | TTree cTree(name,"ITS clusters V2"); | |
240 | TClonesArray *clusters=new TClonesArray("AliITSclusterV2",1000); | |
241 | cTree.Branch("Clusters",&clusters); | |
242 | ||
243 | static TClonesArray *points=its->RecPoints(); | |
244 | AliITSsimulationFastPoints sim; | |
245 | Int_t ncl=0; | |
246 | for (Int_t m=0; m<mmax; m++) { | |
247 | AliITSmodule *mod=its->GetModule(m); | |
248 | sim.CreateFastRecPoints(mod,m,gRandom); | |
249 | ||
250 | Int_t s = (m<=fLastSPD1) ? -1 : 1; | |
251 | RecPoints2Clusters(points, fYshift[m], fZshift[m], fNdet[m], s, clusters); | |
252 | its->ResetRecPoints(); | |
253 | ||
254 | ncl+=clusters->GetEntriesFast(); | |
255 | cTree.Fill(); | |
256 | clusters->Clear(); | |
257 | } | |
258 | cTree.Write(); | |
259 | ||
260 | cerr<<"Number of found fast clusters : "<<ncl<<endl; | |
261 | ||
262 | delete clusters; | |
263 | ||
264 | savedir->cd(); | |
265 | } | |
266 | ||
267 | //*********************************** | |
268 | ||
269 | #ifndef V1 | |
270 | ||
271 | void AliITSclustererV2:: | |
272 | FindCluster(Int_t k,Int_t maxz,AliBin *bins,Int_t &n,Int_t *idx) { | |
273 | //------------------------------------------------------------ | |
274 | // returns an array of indices of digits belonging to the cluster | |
275 | // (needed when the segmentation is not regular) | |
276 | //------------------------------------------------------------ | |
277 | if (n<200) idx[n++]=bins[k].GetIndex(); | |
278 | bins[k].Use(); | |
279 | ||
280 | if (bins[k-maxz].IsNotUsed()) FindCluster(k-maxz,maxz,bins,n,idx); | |
281 | if (bins[k-1 ].IsNotUsed()) FindCluster(k-1 ,maxz,bins,n,idx); | |
282 | if (bins[k+maxz].IsNotUsed()) FindCluster(k+maxz,maxz,bins,n,idx); | |
283 | if (bins[k+1 ].IsNotUsed()) FindCluster(k+1 ,maxz,bins,n,idx); | |
284 | /* | |
285 | if (bins[k-maxz-1].IsNotUsed()) FindCluster(k-maxz-1,maxz,bins,n,idx); | |
286 | if (bins[k-maxz+1].IsNotUsed()) FindCluster(k-maxz+1,maxz,bins,n,idx); | |
287 | if (bins[k+maxz-1].IsNotUsed()) FindCluster(k+maxz-1,maxz,bins,n,idx); | |
288 | if (bins[k+maxz+1].IsNotUsed()) FindCluster(k+maxz+1,maxz,bins,n,idx); | |
289 | */ | |
290 | } | |
291 | ||
292 | void AliITSclustererV2:: | |
293 | FindClustersSPD(const TClonesArray *digits, TClonesArray *clusters) { | |
294 | //------------------------------------------------------------ | |
295 | // Actual SPD cluster finder | |
296 | //------------------------------------------------------------ | |
297 | const Int_t kMAXBIN=(fNzSPD+2)*(fNySPD+2); | |
298 | ||
299 | Int_t ndigits=digits->GetEntriesFast(); | |
300 | AliBin *bins=new AliBin[kMAXBIN]; | |
301 | ||
302 | Int_t k; | |
303 | AliITSdigitSPD *d=0; | |
304 | for (k=0; k<ndigits; k++) { | |
305 | d=(AliITSdigitSPD*)digits->UncheckedAt(k); | |
306 | Int_t i=d->GetCoord2()+1; //y | |
307 | Int_t j=d->GetCoord1()+1; | |
308 | bins[i*fNzSPD+j].SetIndex(k); | |
309 | bins[i*fNzSPD+j].SetMask(1); | |
310 | } | |
311 | ||
312 | Int_t n=0; TClonesArray &cl=*clusters; | |
313 | for (k=0; k<kMAXBIN; k++) { | |
314 | if (!bins[k].IsNotUsed()) continue; | |
315 | Int_t ni=0, idx[200]; | |
316 | FindCluster(k,fNzSPD,bins,ni,idx); | |
317 | if (ni==200) {cerr<<"SPD: Too big cluster !\n"; continue;} | |
318 | ||
319 | Int_t lab[4]; | |
320 | lab[0]=-2; | |
321 | lab[1]=-2; | |
322 | lab[2]=-2; | |
323 | lab[3]=fNdet[fI]; | |
324 | ||
325 | d=(AliITSdigitSPD*)digits->UncheckedAt(idx[0]); | |
326 | Int_t ymin=d->GetCoord2(),ymax=ymin; | |
327 | Int_t zmin=d->GetCoord1(),zmax=zmin; | |
328 | Float_t y=0.,z=0.,q=0.; | |
329 | for (Int_t l=0; l<ni; l++) { | |
330 | d=(AliITSdigitSPD*)digits->UncheckedAt(idx[l]); | |
331 | ||
332 | if (ymin > d->GetCoord2()) ymin=d->GetCoord2(); | |
333 | if (ymax < d->GetCoord2()) ymax=d->GetCoord2(); | |
334 | if (zmin > d->GetCoord1()) zmin=d->GetCoord1(); | |
335 | if (zmax < d->GetCoord1()) zmax=d->GetCoord1(); | |
336 | ||
337 | Int_t lab0=(d->GetTracks())[0]; | |
338 | if (lab0>=0) { | |
339 | if (lab[0]<0) { | |
340 | lab[0]=lab0; | |
341 | } else if (lab[1]<0) { | |
342 | if (lab0!=lab[0]) lab[1]=lab0; | |
343 | } else if (lab[2]<0) { | |
344 | if (lab0!=lab[0]) | |
345 | if (lab0!=lab[1]) lab[2]=lab0; | |
346 | } | |
347 | } | |
348 | Float_t qq=d->GetSignal(); | |
349 | y+=qq*fYSPD[d->GetCoord2()]; z+=qq*fZSPD[d->GetCoord1()]; q+=qq; | |
350 | } | |
351 | y/=q; z/=q; | |
352 | y-=fHwSPD; z-=fHlSPD; | |
353 | ||
354 | Float_t lp[5]; | |
355 | lp[0]=-y-fYshift[fI]; if (fI<=fLastSPD1) lp[0]=-lp[0]; | |
356 | lp[1]= z+fZshift[fI]; | |
357 | lp[2]= fYpitchSPD*fYpitchSPD/12.; | |
358 | lp[3]= fZ1pitchSPD*fZ1pitchSPD/12.; | |
359 | //lp[4]= q; | |
360 | lp[4]= (zmax-zmin+1)*100 + (ymax-ymin+1); | |
361 | ||
362 | //CheckLabels(lab); | |
363 | new (cl[n]) AliITSclusterV2(lab,lp); n++; | |
364 | } | |
365 | ||
366 | delete bins; | |
367 | } | |
368 | ||
369 | Bool_t AliITSclustererV2::IsMaximum(Int_t k,Int_t max,const AliBin *bins) { | |
370 | //------------------------------------------------------------ | |
371 | //is this a local maximum ? | |
372 | //------------------------------------------------------------ | |
373 | UShort_t q=bins[k].GetQ(); | |
374 | if (q==1023) return kFALSE; | |
375 | if (bins[k-max].GetQ() > q) return kFALSE; | |
376 | if (bins[k-1 ].GetQ() > q) return kFALSE; | |
377 | if (bins[k+max].GetQ() > q) return kFALSE; | |
378 | if (bins[k+1 ].GetQ() > q) return kFALSE; | |
379 | if (bins[k-max-1].GetQ() > q) return kFALSE; | |
380 | if (bins[k+max-1].GetQ() > q) return kFALSE; | |
381 | if (bins[k+max+1].GetQ() > q) return kFALSE; | |
382 | if (bins[k-max+1].GetQ() > q) return kFALSE; | |
383 | return kTRUE; | |
384 | } | |
385 | ||
386 | void AliITSclustererV2:: | |
387 | FindPeaks(Int_t k,Int_t max,AliBin *b,Int_t *idx,UInt_t *msk,Int_t& n) { | |
388 | //------------------------------------------------------------ | |
389 | //find local maxima | |
390 | //------------------------------------------------------------ | |
391 | if (n<31) | |
392 | if (IsMaximum(k,max,b)) { | |
393 | idx[n]=k; msk[n]=(2<<n); | |
394 | n++; | |
395 | } | |
396 | b[k].SetMask(0); | |
397 | if (b[k-max].GetMask()&1) FindPeaks(k-max,max,b,idx,msk,n); | |
398 | if (b[k-1 ].GetMask()&1) FindPeaks(k-1 ,max,b,idx,msk,n); | |
399 | if (b[k+max].GetMask()&1) FindPeaks(k+max,max,b,idx,msk,n); | |
400 | if (b[k+1 ].GetMask()&1) FindPeaks(k+1 ,max,b,idx,msk,n); | |
401 | } | |
402 | ||
403 | void AliITSclustererV2:: | |
404 | MarkPeak(Int_t k, Int_t max, AliBin *bins, UInt_t m) { | |
405 | //------------------------------------------------------------ | |
406 | //mark this peak | |
407 | //------------------------------------------------------------ | |
408 | UShort_t q=bins[k].GetQ(); | |
409 | ||
410 | bins[k].SetMask(bins[k].GetMask()|m); | |
411 | ||
412 | if (bins[k-max].GetQ() <= q) | |
413 | if ((bins[k-max].GetMask()&m) == 0) MarkPeak(k-max,max,bins,m); | |
414 | if (bins[k-1 ].GetQ() <= q) | |
415 | if ((bins[k-1 ].GetMask()&m) == 0) MarkPeak(k-1 ,max,bins,m); | |
416 | if (bins[k+max].GetQ() <= q) | |
417 | if ((bins[k+max].GetMask()&m) == 0) MarkPeak(k+max,max,bins,m); | |
418 | if (bins[k+1 ].GetQ() <= q) | |
419 | if ((bins[k+1 ].GetMask()&m) == 0) MarkPeak(k+1 ,max,bins,m); | |
420 | } | |
421 | ||
422 | void AliITSclustererV2:: | |
423 | MakeCluster(Int_t k,Int_t max,AliBin *bins,UInt_t m,AliITSclusterV2 &c) { | |
424 | //------------------------------------------------------------ | |
425 | //make cluster using digits of this peak | |
426 | //------------------------------------------------------------ | |
427 | Float_t q=(Float_t)bins[k].GetQ(); | |
428 | Int_t i=k/max, j=k-i*max; | |
429 | ||
430 | c.SetQ(c.GetQ()+q); | |
431 | c.SetY(c.GetY()+i*q); | |
432 | c.SetZ(c.GetZ()+j*q); | |
433 | c.SetSigmaY2(c.GetSigmaY2()+i*i*q); | |
434 | c.SetSigmaZ2(c.GetSigmaZ2()+j*j*q); | |
435 | ||
436 | bins[k].SetMask(0xFFFFFFFE); | |
437 | ||
438 | if (bins[k-max].GetMask() == m) MakeCluster(k-max,max,bins,m,c); | |
439 | if (bins[k-1 ].GetMask() == m) MakeCluster(k-1 ,max,bins,m,c); | |
440 | if (bins[k+max].GetMask() == m) MakeCluster(k+max,max,bins,m,c); | |
441 | if (bins[k+1 ].GetMask() == m) MakeCluster(k+1 ,max,bins,m,c); | |
442 | } | |
443 | ||
444 | void AliITSclustererV2:: | |
445 | FindClustersSDD(const TClonesArray *digits, TClonesArray *clusters) { | |
446 | //------------------------------------------------------------ | |
447 | // Actual SDD cluster finder | |
448 | //------------------------------------------------------------ | |
449 | const Int_t kMAXBIN=(fNzSDD+2)*(fNySDD+2); | |
450 | ||
451 | AliBin *bins[2]; | |
452 | bins[0]=new AliBin[kMAXBIN]; | |
453 | bins[1]=new AliBin[kMAXBIN]; | |
454 | ||
455 | AliITSdigitSDD *d=0; | |
456 | Int_t i, ndigits=digits->GetEntriesFast(); | |
457 | for (i=0; i<ndigits; i++) { | |
458 | d=(AliITSdigitSDD*)digits->UncheckedAt(i); | |
459 | Int_t y=d->GetCoord2()+1; //y | |
460 | Int_t z=d->GetCoord1()+1; //z | |
461 | Int_t q=d->GetSignal(); | |
462 | if (z <= fNzSDD) { | |
463 | bins[0][y*fNzSDD+z].SetQ(q); | |
464 | bins[0][y*fNzSDD+z].SetMask(1); | |
465 | bins[0][y*fNzSDD+z].SetIndex(i); | |
466 | } else { | |
467 | z-=fNzSDD; | |
468 | bins[1][y*fNzSDD+z].SetQ(q); | |
469 | bins[1][y*fNzSDD+z].SetMask(1); | |
470 | bins[1][y*fNzSDD+z].SetIndex(i); | |
471 | } | |
472 | } | |
473 | ||
474 | Int_t ncl=0; TClonesArray &cl=*clusters; | |
475 | for (Int_t s=0; s<2; s++) | |
476 | for (i=0; i<kMAXBIN; i++) { | |
477 | if (bins[s][i].IsUsed()) continue; | |
478 | Int_t idx[32]; UInt_t msk[32]; Int_t npeaks=0; | |
479 | FindPeaks(i, fNzSDD, bins[s], idx, msk, npeaks); | |
480 | ||
481 | if (npeaks>30) continue; | |
482 | ||
483 | Int_t k,l; | |
484 | for (k=0; k<npeaks-1; k++){//mark adjacent peaks | |
485 | if (idx[k] < 0) continue; //this peak is already removed | |
486 | for (l=k+1; l<npeaks; l++) { | |
487 | if (idx[l] < 0) continue; //this peak is already removed | |
488 | Int_t ki=idx[k]/fNzSDD, kj=idx[k] - ki*fNzSDD; | |
489 | Int_t li=idx[l]/fNzSDD, lj=idx[l] - li*fNzSDD; | |
490 | Int_t di=TMath::Abs(ki - li); | |
491 | Int_t dj=TMath::Abs(kj - lj); | |
492 | if (di>1 || dj>1) continue; | |
493 | if (bins[s][idx[k]].GetQ() > bins[s][idx[l]].GetQ()) { | |
494 | msk[l]=msk[k]; | |
495 | idx[l]*=-1; | |
496 | } else { | |
497 | msk[k]=msk[l]; | |
498 | idx[k]*=-1; | |
499 | break; | |
500 | } | |
501 | } | |
502 | } | |
503 | ||
504 | for (k=0; k<npeaks; k++) { | |
505 | MarkPeak(TMath::Abs(idx[k]), fNzSDD, bins[s], msk[k]); | |
506 | } | |
507 | ||
508 | for (k=0; k<npeaks; k++) { | |
509 | if (idx[k] < 0) continue; //removed peak | |
510 | AliITSclusterV2 c; | |
511 | MakeCluster(idx[k], fNzSDD, bins[s], msk[k], c); | |
512 | ||
513 | //if (c.GetQ() < 200) continue; //noise cluster | |
514 | ||
515 | /* | |
516 | Float_t s2 = c.GetSigmaY2()/c.GetQ() - c.GetY()*c.GetY(); | |
517 | Float_t w=par->GetPadPitchWidth(sec); | |
518 | c.SetSigmaY2((s2 + 1./12.)*w*w); | |
519 | if (s2 != 0.) { | |
520 | c.SetSigmaY2(c.GetSigmaY2()*0.108); | |
521 | if (sec<par->GetNInnerSector()) c.SetSigmaY2(c.GetSigmaY2()*2.07); | |
522 | } | |
523 | ||
524 | s2 = c.GetSigmaZ2()/c.GetQ() - c.GetZ()*c.GetZ(); | |
525 | w=par->GetZWidth(); | |
526 | c.SetSigmaZ2((s2 + 1./12.)*w*w); | |
527 | if (s2 != 0.) { | |
528 | c.SetSigmaZ2(c.GetSigmaZ2()*0.169); | |
529 | if (sec<par->GetNInnerSector()) c.SetSigmaZ2(c.GetSigmaZ2()*1.77); | |
530 | } | |
531 | */ | |
532 | ||
533 | c.SetSigmaY2(0.0030*0.0030); | |
534 | c.SetSigmaZ2(0.0020*0.0020); | |
535 | c.SetDetectorIndex(fNdet[fI]); | |
536 | ||
537 | Float_t y=c.GetY(),z=c.GetZ(), q=c.GetQ(); | |
538 | y/=q; z/=q; | |
539 | ||
540 | y=(y-0.5)*fYpitchSDD; | |
541 | y-=fHwSDD; | |
542 | y-=fYoffSDD; //delay ? | |
543 | if (s) y=-y; | |
544 | ||
545 | z=(z-0.5)*fZpitchSDD; | |
546 | z-=fHlSDD; | |
547 | ||
548 | y=-y-fYshift[fI]; | |
549 | z= z+fZshift[fI]; | |
550 | c.SetY(y); | |
551 | c.SetZ(z); | |
552 | ||
553 | c.SetQ(q/20.); //to be consistent with the SSD charges | |
554 | ||
555 | AliBin *b=&bins[s][idx[k]]; | |
556 | d=(AliITSdigitSDD*)digits->UncheckedAt(b->GetIndex()); | |
557 | Int_t l0=(d->GetTracks())[0]; | |
558 | if (l0<0) { | |
559 | b=&bins[s][idx[k]-1]; | |
560 | if (b->GetQ()>0) { | |
561 | d=(AliITSdigitSDD*)digits->UncheckedAt(b->GetIndex()); | |
562 | l0=(d->GetTracks())[0]; | |
563 | } | |
564 | } | |
565 | if (l0<0) { | |
566 | b=&bins[s][idx[k]+1]; | |
567 | if (b->GetQ()>0) { | |
568 | d=(AliITSdigitSDD*)digits->UncheckedAt(b->GetIndex()); | |
569 | l0=(d->GetTracks())[0]; | |
570 | } | |
571 | } | |
572 | if (l0<0) { | |
573 | b=&bins[s][idx[k]-fNzSDD]; | |
574 | if (b->GetQ()>0) { | |
575 | d=(AliITSdigitSDD*)digits->UncheckedAt(b->GetIndex()); | |
576 | l0=(d->GetTracks())[0]; | |
577 | } | |
578 | } | |
579 | if (l0<0) { | |
580 | b=&bins[s][idx[k]+fNzSDD]; | |
581 | if (b->GetQ()>0) { | |
582 | d=(AliITSdigitSDD*)digits->UncheckedAt(b->GetIndex()); | |
583 | l0=(d->GetTracks())[0]; | |
584 | } | |
585 | } | |
586 | ||
587 | if (l0<0) { | |
588 | b=&bins[s][idx[k]+fNzSDD+1]; | |
589 | if (b->GetQ()>0) { | |
590 | d=(AliITSdigitSDD*)digits->UncheckedAt(b->GetIndex()); | |
591 | l0=(d->GetTracks())[0]; | |
592 | } | |
593 | } | |
594 | if (l0<0) { | |
595 | b=&bins[s][idx[k]+fNzSDD-1]; | |
596 | if (b->GetQ()>0) { | |
597 | d=(AliITSdigitSDD*)digits->UncheckedAt(b->GetIndex()); | |
598 | l0=(d->GetTracks())[0]; | |
599 | } | |
600 | } | |
601 | if (l0<0) { | |
602 | b=&bins[s][idx[k]-fNzSDD+1]; | |
603 | if (b->GetQ()>0) { | |
604 | d=(AliITSdigitSDD*)digits->UncheckedAt(b->GetIndex()); | |
605 | l0=(d->GetTracks())[0]; | |
606 | } | |
607 | } | |
608 | if (l0<0) { | |
609 | b=&bins[s][idx[k]-fNzSDD-1]; | |
610 | if (b->GetQ()>0) { | |
611 | d=(AliITSdigitSDD*)digits->UncheckedAt(b->GetIndex()); | |
612 | l0=(d->GetTracks())[0]; | |
613 | } | |
614 | } | |
615 | ||
616 | { | |
617 | Int_t lab[3]; | |
618 | lab[0]=(d->GetTracks())[0]; | |
619 | lab[1]=(d->GetTracks())[1]; | |
620 | lab[2]=(d->GetTracks())[2]; | |
621 | //CheckLabels(lab); | |
622 | c.SetLabel(lab[0],0); | |
623 | c.SetLabel(lab[1],1); | |
624 | c.SetLabel(lab[2],2); | |
625 | } | |
626 | ||
627 | new (cl[ncl]) AliITSclusterV2(c); ncl++; | |
628 | } | |
629 | } | |
630 | ||
631 | delete[] bins[0]; | |
632 | delete[] bins[1]; | |
633 | } | |
634 | ||
635 | void AliITSclustererV2:: | |
636 | FindClustersSSD(const TClonesArray *digits, TClonesArray *clusters) { | |
637 | //------------------------------------------------------------ | |
638 | // Actual SSD cluster finder | |
639 | //------------------------------------------------------------ | |
640 | Int_t smax=digits->GetEntriesFast(); | |
641 | if (smax==0) return; | |
642 | ||
643 | const Int_t MAX=1000; | |
644 | Int_t np=0, nn=0; | |
645 | Ali1Dcluster pos[MAX], neg[MAX]; | |
646 | Float_t y=0., q=0., qmax=0.; | |
647 | Int_t lab[4]={-2,-2,-2,-2}; | |
648 | ||
649 | TClonesArray &cl=*clusters; | |
650 | ||
651 | AliITSdigitSSD *d=(AliITSdigitSSD*)digits->UncheckedAt(0); | |
652 | q += d->GetSignal(); | |
653 | y += d->GetCoord2()*d->GetSignal(); | |
654 | qmax=d->GetSignal(); | |
655 | lab[0]=d->GetTrack(0); lab[1]=d->GetTrack(1); lab[2]=d->GetTrack(2); | |
656 | Int_t curr=d->GetCoord2(); | |
657 | Int_t flag=d->GetCoord1(); | |
658 | Int_t *n=&nn; | |
659 | Ali1Dcluster *c=neg; | |
660 | Int_t nd=0; | |
661 | for (Int_t s=1; s<smax; s++) { | |
662 | d=(AliITSdigitSSD*)digits->UncheckedAt(s); | |
663 | Int_t strip=d->GetCoord2(); | |
664 | if ((strip-curr) > 1 || flag!=d->GetCoord1()) { | |
665 | c[*n].SetY(y/q); | |
666 | c[*n].SetQ(q); | |
667 | c[*n].SetNd(nd); | |
668 | c[*n].SetLabels(lab); | |
669 | //Split suspiciously big cluster | |
670 | if (nd>3) { | |
671 | c[*n].SetY(y/q-0.5*nd); | |
672 | c[*n].SetQ(0.5*q); | |
673 | (*n)++; | |
674 | if (*n==MAX) { | |
675 | cerr<< | |
676 | "AliITSclustererV2::FindClustersSSD: Too many 1D clusters !\n"; | |
677 | return; | |
678 | } | |
679 | c[*n].SetY(y/q+0.5*nd); | |
680 | c[*n].SetQ(0.5*q); | |
681 | c[*n].SetNd(nd); | |
682 | c[*n].SetLabels(lab); | |
683 | } | |
684 | (*n)++; | |
685 | if (*n==MAX) { | |
686 | cerr<<"AliITSclustererV2::FindClustersSSD: Too many 1D clusters !\n"; | |
687 | return; | |
688 | } | |
689 | y=q=qmax=0.; | |
690 | nd=0; | |
691 | lab[0]=lab[1]=lab[2]=-2; | |
692 | if (flag!=d->GetCoord1()) { n=&np; c=pos; } | |
693 | } | |
694 | flag=d->GetCoord1(); | |
695 | q += d->GetSignal(); | |
696 | y += d->GetCoord2()*d->GetSignal(); | |
697 | nd++; | |
698 | if (d->GetSignal()>qmax) { | |
699 | qmax=d->GetSignal(); | |
700 | lab[0]=d->GetTrack(0); lab[1]=d->GetTrack(1); lab[2]=d->GetTrack(2); | |
701 | } | |
702 | curr=strip; | |
703 | } | |
704 | c[*n].SetY(y/q); | |
705 | c[*n].SetQ(q); | |
706 | c[*n].SetNd(nd); | |
707 | c[*n].SetLabels(lab); | |
708 | //Split suspiciously big cluster | |
709 | if (nd>3) { | |
710 | c[*n].SetY(y/q-0.5*nd); | |
711 | c[*n].SetQ(0.5*q); | |
712 | (*n)++; | |
713 | if (*n==MAX) { | |
714 | cerr<<"AliITSclustererV2::FindClustersSSD: Too many 1D clusters !\n"; | |
715 | return; | |
716 | } | |
717 | c[*n].SetY(y/q+0.5*nd); | |
718 | c[*n].SetQ(0.5*q); | |
719 | c[*n].SetNd(nd); | |
720 | c[*n].SetLabels(lab); | |
721 | } | |
722 | (*n)++; | |
723 | if (*n==MAX) { | |
724 | cerr<<"AliITSclustererV2::FindClustersSSD: Too many 1D clusters !\n"; | |
725 | return; | |
726 | } | |
727 | ||
728 | Float_t tanp=fTanP, tann=fTanN; | |
729 | if (fI>fLastSSD1) {tann=fTanP; tanp=fTanN;} | |
730 | ||
731 | Int_t idet=fNdet[fI]; | |
732 | Int_t ncl=0; | |
733 | for (Int_t i=0; i<np; i++) { | |
734 | //Float_t dq_min=1.e+33; | |
735 | Float_t ybest=1000,zbest=1000,qbest=0; | |
736 | Float_t yp=pos[i].GetY()*fYpitchSSD; | |
737 | for (Int_t j=0; j<nn; j++) { | |
738 | //if (pos[i].fTracks[0] != neg[j].fTracks[0]) continue; | |
739 | Float_t yn=neg[j].GetY()*fYpitchSSD; | |
740 | Float_t zt=(2*fHlSSD*tanp + yp - yn)/(tann+tanp); | |
741 | Float_t yt=yn + tann*zt; | |
742 | zt-=fHlSSD; yt-=fHwSSD; | |
743 | if (TMath::Abs(yt)<fHwSSD+0.01) | |
744 | if (TMath::Abs(zt)<fHlSSD+0.01) { | |
745 | //if (TMath::Abs(pos[i].GetQ()-neg[j].GetQ())<dq_min) { | |
746 | //dq_min=TMath::Abs(pos[i].GetQ()-neg[j].GetQ()); | |
747 | ybest=yt; zbest=zt; | |
748 | qbest=0.5*(pos[i].GetQ()+neg[j].GetQ()); | |
749 | ||
750 | lab[0]=pos[i].GetLabel(0); | |
751 | lab[1]=pos[i].GetLabel(1); | |
752 | lab[2]=neg[i].GetLabel(0); | |
753 | lab[3]=(((i<<10) + j)<<10) + idet; // pos|neg|det | |
754 | Float_t lp[5]; | |
755 | lp[0]=-ybest-fYshift[fI]; | |
756 | lp[1]= zbest+fZshift[fI]; | |
757 | lp[2]=0.0025*0.0025; //SigmaY2 | |
758 | lp[3]=0.110*0.110; //SigmaZ2 | |
759 | if (pos[i].GetNd()+neg[j].GetNd() > 4) { | |
760 | lp[2]*=9; | |
761 | lp[3]*=9; | |
762 | } | |
763 | lp[4]=qbest; //Q | |
764 | ||
765 | //CheckLabels(lab); | |
766 | new (cl[ncl]) AliITSclusterV2(lab,lp); ncl++; | |
767 | } | |
768 | } | |
769 | /* | |
770 | if (ybest<100) { | |
771 | lab[3]=idet; | |
772 | Float_t lp[5]; | |
773 | lp[0]=-ybest-fYshift[fI]; | |
774 | lp[1]= zbest+fZshift[fI]; | |
775 | lp[2]=0.002*0.002; //SigmaY2 | |
776 | lp[3]=0.080*0.080; //SigmaZ2 | |
777 | lp[4]=qbest; //Q | |
778 | // | |
779 | new (cl[ncl]) AliITSclusterV2(lab,lp); ncl++; | |
780 | } | |
781 | */ | |
782 | } | |
783 | } | |
784 | ||
785 | #else //V1 | |
786 | ||
787 | #include "AliITSDetType.h" | |
788 | ||
789 | #include "AliITSsegmentationSPD.h" | |
790 | #include "AliITSClusterFinderSPD.h" | |
791 | ||
792 | #include "AliITSresponseSDD.h" | |
793 | #include "AliITSsegmentationSDD.h" | |
794 | #include "AliITSClusterFinderSDD.h" | |
795 | ||
796 | #include "AliITSsegmentationSSD.h" | |
797 | #include "AliITSClusterFinderSSD.h" | |
798 | ||
799 | ||
800 | void AliITSclustererV2:: | |
801 | FindClustersSPD(const TClonesArray *digits, TClonesArray *clusters) { | |
802 | //------------------------------------------------------------ | |
803 | // Actual SPD cluster finding based on AliITSClusterFinderSPD | |
804 | //------------------------------------------------------------ | |
805 | static AliITS *its=(AliITS*)gAlice->GetModule("ITS"); | |
806 | static TClonesArray *points=its->RecPoints(); | |
807 | static AliITSsegmentationSPD *seg= | |
808 | (AliITSsegmentationSPD *)its->DetType(0)->GetSegmentationModel(); | |
809 | static AliITSClusterFinderSPD cf(seg, (TClonesArray*)digits, points); | |
810 | ||
811 | cf.FindRawClusters(fI); | |
812 | Int_t s = (fI<=fLastSPD1) ? -1 : 1; | |
813 | RecPoints2Clusters(points, fYshift[fI], fZshift[fI], fNdet[fI], s, clusters); | |
814 | its->ResetRecPoints(); | |
815 | ||
816 | } | |
817 | ||
818 | void AliITSclustererV2:: | |
819 | FindClustersSDD(const TClonesArray *digits, TClonesArray *clusters) { | |
820 | //------------------------------------------------------------ | |
821 | // Actual SDD cluster finding based on AliITSClusterFinderSDD | |
822 | //------------------------------------------------------------ | |
823 | static AliITS *its=(AliITS*)gAlice->GetModule("ITS"); | |
824 | static TClonesArray *points=its->RecPoints(); | |
825 | static AliITSresponseSDD *resp= | |
826 | (AliITSresponseSDD *)its->DetType(1)->GetResponseModel(); | |
827 | static AliITSsegmentationSDD *seg= | |
828 | (AliITSsegmentationSDD *)its->DetType(1)->GetSegmentationModel(); | |
829 | static AliITSClusterFinderSDD | |
830 | cf(seg,resp,(TClonesArray*)digits,its->ClustersAddress(1)); | |
831 | ||
832 | cf.FindRawClusters(fI); | |
833 | RecPoints2Clusters(points, fYshift[fI], fZshift[fI], fNdet[fI], 1, clusters); | |
834 | its->ResetClusters(1); | |
835 | its->ResetRecPoints(); | |
836 | ||
837 | } | |
838 | ||
839 | void AliITSclustererV2:: | |
840 | FindClustersSSD(const TClonesArray *digits, TClonesArray *clusters) { | |
841 | //------------------------------------------------------------ | |
842 | // Actual SSD cluster finding based on AliITSClusterFinderSSD | |
843 | //------------------------------------------------------------ | |
844 | static AliITS *its=(AliITS*)gAlice->GetModule("ITS"); | |
845 | static TClonesArray *points=its->RecPoints(); | |
846 | static AliITSsegmentationSSD *seg= | |
847 | (AliITSsegmentationSSD *)its->DetType(2)->GetSegmentationModel(); | |
848 | static AliITSClusterFinderSSD cf(seg,(TClonesArray*)digits); | |
849 | ||
850 | cf.FindRawClusters(fI); | |
851 | RecPoints2Clusters(points, fYshift[fI], fZshift[fI], fNdet[fI], 1, clusters); | |
852 | its->ResetRecPoints(); | |
853 | ||
854 | } | |
855 | ||
856 | #endif |