]>
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 | ||
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 | ||
23 | ClassImp(AliITSclustererV2) | |
24 | ||
25 | extern AliRun *gAlice; | |
26 | ||
27 | AliITSclustererV2::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 | ||
86 | void 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 | ||
167 | static 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 | 190 | void 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 | ||
214 | void 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 | ||
276 | void AliITSclustererV2:: | |
277 | FindCluster(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 | ||
297 | void AliITSclustererV2:: | |
298 | FindClustersSPD(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 | ||
374 | Bool_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 | ||
391 | void AliITSclustererV2:: | |
392 | FindPeaks(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 | ||
408 | void AliITSclustererV2:: | |
409 | MarkPeak(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 | ||
427 | void AliITSclustererV2:: | |
428 | MakeCluster(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 | ||
449 | void AliITSclustererV2:: | |
450 | FindClustersSDD(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 | ||
640 | void AliITSclustererV2:: | |
641 | FindClustersSSD(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 | ||
805 | void AliITSclustererV2:: | |
806 | FindClustersSPD(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 | ||
822 | void AliITSclustererV2:: | |
823 | FindClustersSDD(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 | ||
848 | void AliITSclustererV2:: | |
849 | FindClustersSSD(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 |