]>
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" | |
c391f9d9 | 15 | #include "AliITSRawStreamSPD.h" |
16 | #include "AliITSRawStreamSDD.h" | |
17 | #include "AliITSRawStreamSSD.h" | |
7f4044f0 | 18 | |
19 | #include <Riostream.h> | |
20 | #include <TFile.h> | |
21 | #include <TTree.h> | |
22 | #include <TClonesArray.h> | |
23 | #include "AliITSgeom.h" | |
24 | #include "AliITSdigit.h" | |
25 | ||
26 | ClassImp(AliITSclustererV2) | |
27 | ||
28 | extern AliRun *gAlice; | |
29 | ||
30 | AliITSclustererV2::AliITSclustererV2(const AliITSgeom *geom) { | |
31 | //------------------------------------------------------------ | |
32 | // Standard constructor | |
33 | //------------------------------------------------------------ | |
34 | AliITSgeom *g=(AliITSgeom*)geom; | |
35 | ||
36 | fEvent=0; | |
37 | fI=0; | |
38 | ||
39 | Int_t mmax=geom->GetIndexMax(); | |
40 | if (mmax>2200) {cerr<<"Too many ITS subdetectors !\n"; exit(1);} | |
41 | Int_t m; | |
42 | for (m=0; m<mmax; m++) { | |
43 | Int_t lay,lad,det; g->GetModuleId(m,lay,lad,det); | |
44 | Float_t x,y,z; g->GetTrans(lay,lad,det,x,y,z); | |
45 | Double_t rot[9]; g->GetRotMatrix(lay,lad,det,rot); | |
46 | fYshift[m] = x*rot[0] + y*rot[1]; | |
47 | fZshift[m] = (Double_t)z; | |
48 | fNdet[m] = (lad-1)*g->GetNdetectors(lay) + (det-1); | |
49 | } | |
c391f9d9 | 50 | fNModules = g->GetIndexMax(); |
7f4044f0 | 51 | |
52 | //SPD geometry | |
53 | fLastSPD1=g->GetModuleIndex(2,1,1)-1; | |
54 | fNySPD=256; fNzSPD=160; | |
55 | fYpitchSPD=0.0050; | |
56 | fZ1pitchSPD=0.0425; fZ2pitchSPD=0.0625; | |
57 | fHwSPD=0.64; fHlSPD=3.48; | |
58 | fYSPD[0]=0.5*fYpitchSPD; | |
59 | for (m=1; m<fNySPD; m++) fYSPD[m]=fYSPD[m-1]+fYpitchSPD; | |
60 | fZSPD[0]=fZ1pitchSPD; | |
61 | for (m=1; m<fNzSPD; m++) { | |
62 | Double_t dz=fZ1pitchSPD; | |
63 | if (m==31 || m==32 || m==63 || m==64 || m==95 || m==96 || | |
64 | m==127 || m==128) dz=fZ2pitchSPD; | |
65 | fZSPD[m]=fZSPD[m-1]+dz; | |
66 | } | |
67 | for (m=0; m<fNzSPD; m++) { | |
68 | Double_t dz=0.5*fZ1pitchSPD; | |
69 | if (m==31 || m==32 || m==63 || m==64 || m==95 || m==96 || | |
70 | m==127 || m==128) dz=0.5*fZ2pitchSPD; | |
71 | fZSPD[m]-=dz; | |
72 | } | |
73 | ||
74 | //SDD geometry | |
75 | fNySDD=256; fNzSDD=256; | |
76 | fYpitchSDD=0.01825; | |
77 | fZpitchSDD=0.02940; | |
78 | fHwSDD=3.5085; fHlSDD=3.7632; | |
79 | fYoffSDD=0.0425; | |
80 | ||
81 | //SSD geometry | |
82 | fLastSSD1=g->GetModuleIndex(6,1,1)-1; | |
83 | fYpitchSSD=0.0095; | |
84 | fHwSSD=3.65; | |
85 | fHlSSD=2.00; | |
86 | fTanP=0.0275; | |
87 | fTanN=0.0075; | |
88 | } | |
89 | ||
90 | void AliITSclustererV2::Digits2Clusters(const TFile *in, TFile *out) { | |
91 | //------------------------------------------------------------ | |
92 | // This function creates ITS clusters | |
93 | //------------------------------------------------------------ | |
94 | Int_t ncl=0; | |
95 | TDirectory *savedir=gDirectory; | |
96 | ||
97 | if (!out->IsOpen()) { | |
98 | cerr<<"AliITSclustererV2::Digits2Clusters(): output file not open !\n"; | |
99 | return; | |
100 | } | |
101 | ||
102 | Char_t name[100]; | |
103 | sprintf(name,"TreeD%d",fEvent); | |
104 | ||
105 | //TTree *dTree=(TTree*)((TFile*)in)->Get(name); | |
106 | TTree *dTree=gAlice->TreeD(); | |
107 | if (!dTree) { | |
108 | cerr<<"Input tree "<<name<<" not found !\n"; | |
109 | return; | |
110 | } | |
111 | ||
112 | TClonesArray *digitsSPD=new TClonesArray("AliITSdigitSPD",3000); | |
113 | dTree->SetBranchAddress("ITSDigitsSPD",&digitsSPD); | |
114 | TClonesArray *digitsSDD=new TClonesArray("AliITSdigitSDD",3000); | |
115 | dTree->SetBranchAddress("ITSDigitsSDD",&digitsSDD); | |
116 | TClonesArray *digitsSSD=new TClonesArray("AliITSdigitSSD",3000); | |
117 | dTree->SetBranchAddress("ITSDigitsSSD",&digitsSSD); | |
118 | ||
119 | Int_t mmax=(Int_t)dTree->GetEntries(); | |
120 | ||
121 | out->cd(); | |
122 | ||
123 | sprintf(name,"TreeC_ITS_%d",fEvent); | |
124 | TTree cTree(name,"ITS clusters V2"); | |
125 | TClonesArray *clusters=new TClonesArray("AliITSclusterV2",1000); | |
126 | cTree.Branch("Clusters",&clusters); | |
127 | ||
128 | for (fI=0; fI<mmax; fI++) { | |
129 | dTree->GetEvent(fI); | |
130 | ||
131 | if (digitsSPD->GetEntriesFast()!=0) | |
132 | FindClustersSPD(digitsSPD,clusters); | |
133 | else if(digitsSDD->GetEntriesFast()!=0) | |
134 | FindClustersSDD(digitsSDD,clusters); | |
135 | else if(digitsSSD->GetEntriesFast()!=0) | |
136 | FindClustersSSD(digitsSSD,clusters); | |
137 | ||
138 | ncl+=clusters->GetEntriesFast(); | |
139 | ||
140 | cTree.Fill(); | |
141 | ||
142 | digitsSPD->Clear(); | |
143 | digitsSDD->Clear(); | |
144 | digitsSSD->Clear(); | |
145 | clusters->Clear(); | |
146 | } | |
147 | cTree.Write(); | |
148 | ||
149 | delete clusters; | |
150 | ||
151 | delete digitsSPD; | |
152 | delete digitsSDD; | |
153 | delete digitsSSD; | |
154 | ||
155 | //delete dTree; | |
156 | ||
157 | cerr<<"Number of found clusters : "<<ncl<<endl; | |
158 | ||
159 | savedir->cd(); | |
160 | } | |
161 | ||
c391f9d9 | 162 | void AliITSclustererV2::Digits2Clusters(TFile *out) { |
163 | //------------------------------------------------------------ | |
164 | // This function creates ITS clusters from raw data | |
165 | //------------------------------------------------------------ | |
166 | TDirectory *savedir=gDirectory; | |
167 | if (!out->IsOpen()) { | |
168 | cerr<<"AliITSclustererV2::Digits2Clusters(): output file not open !\n"; | |
169 | return; | |
170 | } | |
171 | out->cd(); | |
172 | ||
173 | Char_t name[100]; | |
174 | sprintf(name,"TreeC_ITS_%d",fEvent); | |
175 | TTree cTree(name,"ITS clusters V2"); | |
176 | TClonesArray *array=new TClonesArray("AliITSclusterV2",1000); | |
177 | cTree.Branch("Clusters",&array); | |
178 | delete array; | |
179 | ||
f1e0c1c5 | 180 | TClonesArray** clusters = new TClonesArray*[fNModules]; |
c391f9d9 | 181 | // one TClonesArray per module |
182 | ||
183 | AliITSRawStreamSPD inputSPD; | |
184 | FindClustersSPD(&inputSPD, clusters); | |
185 | AliITSRawStreamSDD inputSDD; | |
186 | FindClustersSDD(&inputSDD, clusters); | |
187 | AliITSRawStreamSSD inputSSD; | |
188 | FindClustersSSD(&inputSSD, clusters); | |
189 | ||
190 | // write all clusters to the tree | |
191 | Int_t nClusters = 0; | |
192 | for (Int_t iModule = 0; iModule < fNModules; iModule++) { | |
193 | TClonesArray* array = clusters[iModule]; | |
194 | if (!array) { | |
195 | Error("Digits2Clusters", "data for module %d missing!", iModule); | |
196 | array = new TClonesArray("AliITSclusterV2"); | |
197 | } | |
198 | cTree.SetBranchAddress("Clusters", &array); | |
199 | cTree.Fill(); | |
200 | nClusters += array->GetEntriesFast(); | |
201 | delete array; | |
202 | } | |
203 | cTree.Write(); | |
204 | ||
205 | savedir->cd(); | |
206 | ||
207 | Info("Digits2Clusters", "total number of found clusters in ITS: %d\n", | |
208 | nClusters); | |
209 | } | |
210 | ||
7f4044f0 | 211 | //**** Fast clusters ******************************* |
212 | #include "TParticle.h" | |
213 | ||
214 | #include "AliITS.h" | |
215 | #include "AliITSmodule.h" | |
216 | #include "AliITSRecPoint.h" | |
217 | #include "AliITSsimulationFastPoints.h" | |
218 | #include "AliITSRecPoint.h" | |
219 | ||
220 | static void CheckLabels(Int_t lab[3]) { | |
221 | //------------------------------------------------------------ | |
222 | // Tries to find mother's labels | |
223 | //------------------------------------------------------------ | |
224 | Int_t label=lab[0]; | |
225 | if (label>=0) { | |
226 | TParticle *part=(TParticle*)gAlice->Particle(label); | |
227 | label=-3; | |
228 | while (part->P() < 0.005) { | |
229 | Int_t m=part->GetFirstMother(); | |
230 | if (m<0) {cerr<<"Primary momentum: "<<part->P()<<endl; break;} | |
1cca57bf | 231 | if (part->GetStatusCode()>0) { |
232 | cerr<<"Primary momentum: "<<part->P()<<endl; break; | |
233 | } | |
7f4044f0 | 234 | label=m; |
235 | part=(TParticle*)gAlice->Particle(label); | |
236 | } | |
237 | if (lab[1]<0) lab[1]=label; | |
238 | else if (lab[2]<0) lab[2]=label; | |
239 | else ;//cerr<<"CheckLabels : No empty labels !\n"; | |
240 | } | |
241 | } | |
242 | ||
1cca57bf | 243 | void AliITSclustererV2::RecPoints2Clusters |
244 | (const TClonesArray *points, Int_t idx, TClonesArray *clusters) { | |
7f4044f0 | 245 | //------------------------------------------------------------ |
1cca57bf | 246 | // Conversion AliITSRecPoint -> AliITSclusterV2 for the ITS |
247 | // subdetector indexed by idx | |
7f4044f0 | 248 | //------------------------------------------------------------ |
249 | TClonesArray &cl=*clusters; | |
250 | Int_t ncl=points->GetEntriesFast(); | |
251 | for (Int_t i=0; i<ncl; i++) { | |
252 | AliITSRecPoint *p = (AliITSRecPoint *)points->UncheckedAt(i); | |
253 | Float_t lp[5]; | |
1cca57bf | 254 | lp[0]=-p->GetX()-fYshift[idx]; if (idx<=fLastSPD1) lp[0]*=-1; //SPD1 |
255 | lp[1]=p->GetZ()+fZshift[idx]; | |
256 | lp[2]=p->GetSigmaX2(); | |
257 | lp[3]=p->GetSigmaZ2(); | |
258 | lp[4]=p->GetQ(); | |
7f4044f0 | 259 | Int_t lab[4]; |
260 | lab[0]=p->GetLabel(0); lab[1]=p->GetLabel(1); lab[2]=p->GetLabel(2); | |
1cca57bf | 261 | lab[3]=fNdet[idx]; |
7f4044f0 | 262 | CheckLabels(lab); |
263 | new (cl[i]) AliITSclusterV2(lab,lp); | |
264 | } | |
265 | } | |
266 | ||
267 | void AliITSclustererV2::Hits2Clusters(const TFile *in, TFile *out) { | |
268 | //------------------------------------------------------------ | |
269 | // This function creates ITS clusters | |
270 | //------------------------------------------------------------ | |
271 | TDirectory *savedir=gDirectory; | |
272 | ||
273 | if (!out->IsOpen()) { | |
274 | cerr<<"AliITSclustererV2::Hits2Clusters: output file not open !\n"; | |
275 | return; | |
276 | } | |
277 | ||
278 | if (!gAlice) { | |
279 | cerr<<"AliITSclustererV2::Hits2Clusters : gAlice==0 !\n"; | |
280 | return; | |
281 | } | |
282 | ||
283 | AliITS *its = (AliITS*)gAlice->GetModule("ITS"); | |
284 | if (!its) { | |
285 | cerr<<"AliITSclustererV2::Hits2Clusters : Can't find the ITS !\n"; | |
286 | return; | |
287 | } | |
288 | AliITSgeom *geom=its->GetITSgeom(); | |
289 | Int_t mmax=geom->GetIndexMax(); | |
290 | ||
291 | its->InitModules(-1,mmax); | |
292 | its->FillModules(gAlice->TreeH(),0); | |
293 | ||
294 | out->cd(); | |
295 | ||
296 | Char_t name[100]; | |
297 | sprintf(name,"TreeC_ITS_%d",fEvent); | |
298 | TTree cTree(name,"ITS clusters V2"); | |
299 | TClonesArray *clusters=new TClonesArray("AliITSclusterV2",1000); | |
300 | cTree.Branch("Clusters",&clusters); | |
301 | ||
302 | static TClonesArray *points=its->RecPoints(); | |
303 | AliITSsimulationFastPoints sim; | |
304 | Int_t ncl=0; | |
305 | for (Int_t m=0; m<mmax; m++) { | |
306 | AliITSmodule *mod=its->GetModule(m); | |
307 | sim.CreateFastRecPoints(mod,m,gRandom); | |
308 | ||
1cca57bf | 309 | RecPoints2Clusters(points, m, clusters); |
7f4044f0 | 310 | its->ResetRecPoints(); |
311 | ||
312 | ncl+=clusters->GetEntriesFast(); | |
313 | cTree.Fill(); | |
314 | clusters->Clear(); | |
315 | } | |
316 | cTree.Write(); | |
317 | ||
318 | cerr<<"Number of found fast clusters : "<<ncl<<endl; | |
319 | ||
320 | delete clusters; | |
321 | ||
322 | savedir->cd(); | |
323 | } | |
324 | ||
325 | //*********************************** | |
326 | ||
327 | #ifndef V1 | |
328 | ||
329 | void AliITSclustererV2:: | |
330 | FindCluster(Int_t k,Int_t maxz,AliBin *bins,Int_t &n,Int_t *idx) { | |
331 | //------------------------------------------------------------ | |
332 | // returns an array of indices of digits belonging to the cluster | |
333 | // (needed when the segmentation is not regular) | |
334 | //------------------------------------------------------------ | |
335 | if (n<200) idx[n++]=bins[k].GetIndex(); | |
336 | bins[k].Use(); | |
337 | ||
338 | if (bins[k-maxz].IsNotUsed()) FindCluster(k-maxz,maxz,bins,n,idx); | |
339 | if (bins[k-1 ].IsNotUsed()) FindCluster(k-1 ,maxz,bins,n,idx); | |
340 | if (bins[k+maxz].IsNotUsed()) FindCluster(k+maxz,maxz,bins,n,idx); | |
341 | if (bins[k+1 ].IsNotUsed()) FindCluster(k+1 ,maxz,bins,n,idx); | |
342 | /* | |
343 | if (bins[k-maxz-1].IsNotUsed()) FindCluster(k-maxz-1,maxz,bins,n,idx); | |
344 | if (bins[k-maxz+1].IsNotUsed()) FindCluster(k-maxz+1,maxz,bins,n,idx); | |
345 | if (bins[k+maxz-1].IsNotUsed()) FindCluster(k+maxz-1,maxz,bins,n,idx); | |
346 | if (bins[k+maxz+1].IsNotUsed()) FindCluster(k+maxz+1,maxz,bins,n,idx); | |
347 | */ | |
348 | } | |
349 | ||
350 | void AliITSclustererV2:: | |
351 | FindClustersSPD(const TClonesArray *digits, TClonesArray *clusters) { | |
352 | //------------------------------------------------------------ | |
353 | // Actual SPD cluster finder | |
354 | //------------------------------------------------------------ | |
c391f9d9 | 355 | Int_t kNzBins = fNzSPD + 2; |
356 | const Int_t kMAXBIN=kNzBins*(fNySPD+2); | |
7f4044f0 | 357 | |
358 | Int_t ndigits=digits->GetEntriesFast(); | |
359 | AliBin *bins=new AliBin[kMAXBIN]; | |
360 | ||
361 | Int_t k; | |
362 | AliITSdigitSPD *d=0; | |
363 | for (k=0; k<ndigits; k++) { | |
364 | d=(AliITSdigitSPD*)digits->UncheckedAt(k); | |
365 | Int_t i=d->GetCoord2()+1; //y | |
366 | Int_t j=d->GetCoord1()+1; | |
c391f9d9 | 367 | bins[i*kNzBins+j].SetIndex(k); |
368 | bins[i*kNzBins+j].SetMask(1); | |
7f4044f0 | 369 | } |
370 | ||
371 | Int_t n=0; TClonesArray &cl=*clusters; | |
372 | for (k=0; k<kMAXBIN; k++) { | |
373 | if (!bins[k].IsNotUsed()) continue; | |
374 | Int_t ni=0, idx[200]; | |
c391f9d9 | 375 | FindCluster(k,kNzBins,bins,ni,idx); |
7f4044f0 | 376 | if (ni==200) {cerr<<"SPD: Too big cluster !\n"; continue;} |
377 | ||
378 | Int_t lab[4]; | |
379 | lab[0]=-2; | |
380 | lab[1]=-2; | |
381 | lab[2]=-2; | |
382 | lab[3]=fNdet[fI]; | |
383 | ||
384 | d=(AliITSdigitSPD*)digits->UncheckedAt(idx[0]); | |
385 | Int_t ymin=d->GetCoord2(),ymax=ymin; | |
386 | Int_t zmin=d->GetCoord1(),zmax=zmin; | |
387 | Float_t y=0.,z=0.,q=0.; | |
388 | for (Int_t l=0; l<ni; l++) { | |
389 | d=(AliITSdigitSPD*)digits->UncheckedAt(idx[l]); | |
390 | ||
391 | if (ymin > d->GetCoord2()) ymin=d->GetCoord2(); | |
392 | if (ymax < d->GetCoord2()) ymax=d->GetCoord2(); | |
393 | if (zmin > d->GetCoord1()) zmin=d->GetCoord1(); | |
394 | if (zmax < d->GetCoord1()) zmax=d->GetCoord1(); | |
395 | ||
396 | Int_t lab0=(d->GetTracks())[0]; | |
397 | if (lab0>=0) { | |
398 | if (lab[0]<0) { | |
399 | lab[0]=lab0; | |
400 | } else if (lab[1]<0) { | |
401 | if (lab0!=lab[0]) lab[1]=lab0; | |
402 | } else if (lab[2]<0) { | |
403 | if (lab0!=lab[0]) | |
404 | if (lab0!=lab[1]) lab[2]=lab0; | |
405 | } | |
406 | } | |
407 | Float_t qq=d->GetSignal(); | |
408 | y+=qq*fYSPD[d->GetCoord2()]; z+=qq*fZSPD[d->GetCoord1()]; q+=qq; | |
409 | } | |
410 | y/=q; z/=q; | |
411 | y-=fHwSPD; z-=fHlSPD; | |
412 | ||
413 | Float_t lp[5]; | |
414 | lp[0]=-y-fYshift[fI]; if (fI<=fLastSPD1) lp[0]=-lp[0]; | |
415 | lp[1]= z+fZshift[fI]; | |
416 | lp[2]= fYpitchSPD*fYpitchSPD/12.; | |
417 | lp[3]= fZ1pitchSPD*fZ1pitchSPD/12.; | |
418 | //lp[4]= q; | |
419 | lp[4]= (zmax-zmin+1)*100 + (ymax-ymin+1); | |
420 | ||
421 | //CheckLabels(lab); | |
c391f9d9 | 422 | d=(AliITSdigitSPD*)digits->UncheckedAt(idx[0]); |
7f4044f0 | 423 | new (cl[n]) AliITSclusterV2(lab,lp); n++; |
424 | } | |
425 | ||
426 | delete bins; | |
427 | } | |
428 | ||
c391f9d9 | 429 | void AliITSclustererV2::FindClustersSPD(AliITSRawStream* input, |
430 | TClonesArray** clusters) | |
431 | { | |
432 | //------------------------------------------------------------ | |
433 | // Actual SPD cluster finder for raw data | |
434 | //------------------------------------------------------------ | |
435 | ||
436 | Int_t nClustersSPD = 0; | |
437 | Int_t kNzBins = fNzSPD + 2; | |
438 | Int_t kNyBins = fNySPD + 2; | |
439 | Int_t kMaxBin = kNzBins * kNyBins; | |
440 | AliBin* bins = NULL; | |
441 | ||
442 | // read raw data input stream | |
443 | while (kTRUE) { | |
444 | Bool_t next = input->Next(); | |
445 | if (!next || input->IsNewModule()) { | |
446 | Int_t iModule = input->GetPrevModuleID(); | |
447 | ||
448 | // when all data from a module was read, search for clusters | |
449 | if (bins) { | |
450 | clusters[iModule] = new TClonesArray("AliITSclusterV2"); | |
451 | Int_t nClusters = 0; | |
452 | ||
453 | for (Int_t iBin = 0; iBin < kMaxBin; iBin++) { | |
454 | if (bins[iBin].IsUsed()) continue; | |
455 | Int_t nBins = 0; | |
456 | Int_t idxBins[200]; | |
457 | FindCluster(iBin, kNzBins, bins, nBins, idxBins); | |
458 | if (nBins == 200) { | |
459 | Error("FindClustersSPD", "SPD: Too big cluster !\n"); | |
460 | continue; | |
461 | } | |
462 | ||
463 | Int_t label[4]; | |
464 | label[0] = -2; | |
465 | label[1] = -2; | |
466 | label[2] = -2; | |
467 | // label[3] = iModule; | |
468 | label[3] = fNdet[iModule]; | |
469 | ||
470 | Int_t ymin = (idxBins[0] / kNzBins) - 1; | |
471 | Int_t ymax = ymin; | |
472 | Int_t zmin = (idxBins[0] % kNzBins) - 1; | |
473 | Int_t zmax = zmin; | |
474 | Float_t y = 0.; | |
475 | Float_t z = 0.; | |
476 | Float_t q = 0.; | |
477 | for (Int_t idx = 0; idx < nBins; idx++) { | |
478 | Int_t iy = (idxBins[idx] / kNzBins) - 1; | |
479 | Int_t iz = (idxBins[idx] % kNzBins) - 1; | |
480 | if (ymin > iy) ymin = iy; | |
481 | if (ymax < iy) ymax = iy; | |
482 | if (zmin > iz) zmin = iz; | |
483 | if (zmax < iz) zmax = iz; | |
484 | ||
485 | Float_t qBin = bins[idxBins[idx]].GetQ(); | |
486 | y += qBin * fYSPD[iy]; | |
487 | z += qBin * fZSPD[iz]; | |
488 | q += qBin; | |
489 | } | |
490 | y /= q; | |
491 | z /= q; | |
492 | y -= fHwSPD; | |
493 | z -= fHlSPD; | |
494 | ||
495 | Float_t hit[5]; // y, z, sigma(y)^2, sigma(z)^2, charge | |
496 | hit[0] = -y-fYshift[iModule]; | |
497 | if (iModule <= fLastSPD1) hit[0] = -hit[0]; | |
498 | hit[1] = z+fZshift[iModule]; | |
499 | hit[2] = fYpitchSPD*fYpitchSPD/12.; | |
500 | hit[3] = fZ1pitchSPD*fZ1pitchSPD/12.; | |
501 | // hit[4] = q; | |
502 | hit[4] = (zmax-zmin+1)*100 + (ymax-ymin+1); | |
503 | ||
504 | //CheckLabels(label); | |
505 | new (clusters[iModule]->AddrAt(nClusters)) | |
506 | AliITSclusterV2(label, hit); | |
507 | nClusters++; | |
508 | } | |
509 | ||
510 | nClustersSPD += nClusters; | |
511 | delete bins; | |
512 | } | |
513 | ||
514 | if (!next) break; | |
515 | bins = new AliBin[kMaxBin]; | |
516 | } | |
517 | ||
518 | // fill the current digit into the bins array | |
519 | Int_t index = (input->GetCoord2()+1) * kNzBins + (input->GetCoord1()+1); | |
520 | bins[index].SetIndex(index); | |
521 | bins[index].SetMask(1); | |
522 | bins[index].SetQ(1); | |
523 | } | |
524 | ||
525 | Info("FindClustersSPD", "found clusters in ITS SPD: %d", nClustersSPD); | |
526 | } | |
527 | ||
528 | ||
7f4044f0 | 529 | Bool_t AliITSclustererV2::IsMaximum(Int_t k,Int_t max,const AliBin *bins) { |
530 | //------------------------------------------------------------ | |
531 | //is this a local maximum ? | |
532 | //------------------------------------------------------------ | |
533 | UShort_t q=bins[k].GetQ(); | |
534 | if (q==1023) return kFALSE; | |
535 | if (bins[k-max].GetQ() > q) return kFALSE; | |
536 | if (bins[k-1 ].GetQ() > q) return kFALSE; | |
537 | if (bins[k+max].GetQ() > q) return kFALSE; | |
538 | if (bins[k+1 ].GetQ() > q) return kFALSE; | |
539 | if (bins[k-max-1].GetQ() > q) return kFALSE; | |
540 | if (bins[k+max-1].GetQ() > q) return kFALSE; | |
541 | if (bins[k+max+1].GetQ() > q) return kFALSE; | |
542 | if (bins[k-max+1].GetQ() > q) return kFALSE; | |
543 | return kTRUE; | |
544 | } | |
545 | ||
546 | void AliITSclustererV2:: | |
547 | FindPeaks(Int_t k,Int_t max,AliBin *b,Int_t *idx,UInt_t *msk,Int_t& n) { | |
548 | //------------------------------------------------------------ | |
549 | //find local maxima | |
550 | //------------------------------------------------------------ | |
551 | if (n<31) | |
552 | if (IsMaximum(k,max,b)) { | |
553 | idx[n]=k; msk[n]=(2<<n); | |
554 | n++; | |
555 | } | |
556 | b[k].SetMask(0); | |
557 | if (b[k-max].GetMask()&1) FindPeaks(k-max,max,b,idx,msk,n); | |
558 | if (b[k-1 ].GetMask()&1) FindPeaks(k-1 ,max,b,idx,msk,n); | |
559 | if (b[k+max].GetMask()&1) FindPeaks(k+max,max,b,idx,msk,n); | |
560 | if (b[k+1 ].GetMask()&1) FindPeaks(k+1 ,max,b,idx,msk,n); | |
561 | } | |
562 | ||
563 | void AliITSclustererV2:: | |
564 | MarkPeak(Int_t k, Int_t max, AliBin *bins, UInt_t m) { | |
565 | //------------------------------------------------------------ | |
566 | //mark this peak | |
567 | //------------------------------------------------------------ | |
568 | UShort_t q=bins[k].GetQ(); | |
569 | ||
570 | bins[k].SetMask(bins[k].GetMask()|m); | |
571 | ||
572 | if (bins[k-max].GetQ() <= q) | |
573 | if ((bins[k-max].GetMask()&m) == 0) MarkPeak(k-max,max,bins,m); | |
574 | if (bins[k-1 ].GetQ() <= q) | |
575 | if ((bins[k-1 ].GetMask()&m) == 0) MarkPeak(k-1 ,max,bins,m); | |
576 | if (bins[k+max].GetQ() <= q) | |
577 | if ((bins[k+max].GetMask()&m) == 0) MarkPeak(k+max,max,bins,m); | |
578 | if (bins[k+1 ].GetQ() <= q) | |
579 | if ((bins[k+1 ].GetMask()&m) == 0) MarkPeak(k+1 ,max,bins,m); | |
580 | } | |
581 | ||
582 | void AliITSclustererV2:: | |
583 | MakeCluster(Int_t k,Int_t max,AliBin *bins,UInt_t m,AliITSclusterV2 &c) { | |
584 | //------------------------------------------------------------ | |
585 | //make cluster using digits of this peak | |
586 | //------------------------------------------------------------ | |
587 | Float_t q=(Float_t)bins[k].GetQ(); | |
588 | Int_t i=k/max, j=k-i*max; | |
589 | ||
590 | c.SetQ(c.GetQ()+q); | |
591 | c.SetY(c.GetY()+i*q); | |
592 | c.SetZ(c.GetZ()+j*q); | |
593 | c.SetSigmaY2(c.GetSigmaY2()+i*i*q); | |
594 | c.SetSigmaZ2(c.GetSigmaZ2()+j*j*q); | |
595 | ||
596 | bins[k].SetMask(0xFFFFFFFE); | |
597 | ||
598 | if (bins[k-max].GetMask() == m) MakeCluster(k-max,max,bins,m,c); | |
599 | if (bins[k-1 ].GetMask() == m) MakeCluster(k-1 ,max,bins,m,c); | |
600 | if (bins[k+max].GetMask() == m) MakeCluster(k+max,max,bins,m,c); | |
601 | if (bins[k+1 ].GetMask() == m) MakeCluster(k+1 ,max,bins,m,c); | |
602 | } | |
603 | ||
604 | void AliITSclustererV2:: | |
c391f9d9 | 605 | FindClustersSDD(AliBin* bins[2], Int_t nMaxBin, Int_t nzBins, |
606 | const TClonesArray *digits, TClonesArray *clusters) { | |
7f4044f0 | 607 | //------------------------------------------------------------ |
608 | // Actual SDD cluster finder | |
609 | //------------------------------------------------------------ | |
7f4044f0 | 610 | Int_t ncl=0; TClonesArray &cl=*clusters; |
611 | for (Int_t s=0; s<2; s++) | |
c391f9d9 | 612 | for (Int_t i=0; i<nMaxBin; i++) { |
7f4044f0 | 613 | if (bins[s][i].IsUsed()) continue; |
614 | Int_t idx[32]; UInt_t msk[32]; Int_t npeaks=0; | |
c391f9d9 | 615 | FindPeaks(i, nzBins, bins[s], idx, msk, npeaks); |
7f4044f0 | 616 | |
617 | if (npeaks>30) continue; | |
618 | ||
619 | Int_t k,l; | |
620 | for (k=0; k<npeaks-1; k++){//mark adjacent peaks | |
621 | if (idx[k] < 0) continue; //this peak is already removed | |
622 | for (l=k+1; l<npeaks; l++) { | |
623 | if (idx[l] < 0) continue; //this peak is already removed | |
c391f9d9 | 624 | Int_t ki=idx[k]/nzBins, kj=idx[k] - ki*nzBins; |
625 | Int_t li=idx[l]/nzBins, lj=idx[l] - li*nzBins; | |
7f4044f0 | 626 | Int_t di=TMath::Abs(ki - li); |
627 | Int_t dj=TMath::Abs(kj - lj); | |
628 | if (di>1 || dj>1) continue; | |
629 | if (bins[s][idx[k]].GetQ() > bins[s][idx[l]].GetQ()) { | |
630 | msk[l]=msk[k]; | |
631 | idx[l]*=-1; | |
632 | } else { | |
633 | msk[k]=msk[l]; | |
634 | idx[k]*=-1; | |
635 | break; | |
636 | } | |
637 | } | |
638 | } | |
639 | ||
640 | for (k=0; k<npeaks; k++) { | |
c391f9d9 | 641 | MarkPeak(TMath::Abs(idx[k]), nzBins, bins[s], msk[k]); |
7f4044f0 | 642 | } |
643 | ||
644 | for (k=0; k<npeaks; k++) { | |
645 | if (idx[k] < 0) continue; //removed peak | |
646 | AliITSclusterV2 c; | |
c391f9d9 | 647 | MakeCluster(idx[k], nzBins, bins[s], msk[k], c); |
7f4044f0 | 648 | |
649 | //if (c.GetQ() < 200) continue; //noise cluster | |
650 | ||
651 | /* | |
652 | Float_t s2 = c.GetSigmaY2()/c.GetQ() - c.GetY()*c.GetY(); | |
653 | Float_t w=par->GetPadPitchWidth(sec); | |
654 | c.SetSigmaY2((s2 + 1./12.)*w*w); | |
655 | if (s2 != 0.) { | |
656 | c.SetSigmaY2(c.GetSigmaY2()*0.108); | |
657 | if (sec<par->GetNInnerSector()) c.SetSigmaY2(c.GetSigmaY2()*2.07); | |
658 | } | |
659 | ||
660 | s2 = c.GetSigmaZ2()/c.GetQ() - c.GetZ()*c.GetZ(); | |
661 | w=par->GetZWidth(); | |
662 | c.SetSigmaZ2((s2 + 1./12.)*w*w); | |
663 | if (s2 != 0.) { | |
664 | c.SetSigmaZ2(c.GetSigmaZ2()*0.169); | |
665 | if (sec<par->GetNInnerSector()) c.SetSigmaZ2(c.GetSigmaZ2()*1.77); | |
666 | } | |
667 | */ | |
668 | ||
669 | c.SetSigmaY2(0.0030*0.0030); | |
670 | c.SetSigmaZ2(0.0020*0.0020); | |
671 | c.SetDetectorIndex(fNdet[fI]); | |
672 | ||
673 | Float_t y=c.GetY(),z=c.GetZ(), q=c.GetQ(); | |
674 | y/=q; z/=q; | |
675 | ||
676 | y=(y-0.5)*fYpitchSDD; | |
677 | y-=fHwSDD; | |
678 | y-=fYoffSDD; //delay ? | |
679 | if (s) y=-y; | |
680 | ||
681 | z=(z-0.5)*fZpitchSDD; | |
682 | z-=fHlSDD; | |
683 | ||
684 | y=-y-fYshift[fI]; | |
685 | z= z+fZshift[fI]; | |
686 | c.SetY(y); | |
687 | c.SetZ(z); | |
688 | ||
689 | c.SetQ(q/20.); //to be consistent with the SSD charges | |
690 | ||
c391f9d9 | 691 | if (digits) { |
692 | AliBin *b=&bins[s][idx[k]]; | |
693 | AliITSdigitSDD* d=(AliITSdigitSDD*)digits->UncheckedAt(b->GetIndex()); | |
694 | Int_t l0=(d->GetTracks())[0]; | |
695 | if (l0<0) { | |
696 | b=&bins[s][idx[k]-1]; | |
697 | if (b->GetQ()>0) { | |
698 | d=(AliITSdigitSDD*)digits->UncheckedAt(b->GetIndex()); | |
699 | l0=(d->GetTracks())[0]; | |
700 | } | |
701 | } | |
702 | if (l0<0) { | |
703 | b=&bins[s][idx[k]+1]; | |
704 | if (b->GetQ()>0) { | |
705 | d=(AliITSdigitSDD*)digits->UncheckedAt(b->GetIndex()); | |
706 | l0=(d->GetTracks())[0]; | |
707 | } | |
708 | } | |
709 | if (l0<0) { | |
710 | b=&bins[s][idx[k]-nzBins]; | |
711 | if (b->GetQ()>0) { | |
712 | d=(AliITSdigitSDD*)digits->UncheckedAt(b->GetIndex()); | |
713 | l0=(d->GetTracks())[0]; | |
714 | } | |
715 | } | |
716 | if (l0<0) { | |
717 | b=&bins[s][idx[k]+nzBins]; | |
718 | if (b->GetQ()>0) { | |
719 | d=(AliITSdigitSDD*)digits->UncheckedAt(b->GetIndex()); | |
720 | l0=(d->GetTracks())[0]; | |
721 | } | |
7f4044f0 | 722 | } |
7f4044f0 | 723 | |
c391f9d9 | 724 | if (l0<0) { |
725 | b=&bins[s][idx[k]+nzBins+1]; | |
726 | if (b->GetQ()>0) { | |
727 | d=(AliITSdigitSDD*)digits->UncheckedAt(b->GetIndex()); | |
728 | l0=(d->GetTracks())[0]; | |
729 | } | |
730 | } | |
731 | if (l0<0) { | |
732 | b=&bins[s][idx[k]+nzBins-1]; | |
733 | if (b->GetQ()>0) { | |
734 | d=(AliITSdigitSDD*)digits->UncheckedAt(b->GetIndex()); | |
735 | l0=(d->GetTracks())[0]; | |
736 | } | |
737 | } | |
738 | if (l0<0) { | |
739 | b=&bins[s][idx[k]-nzBins+1]; | |
740 | if (b->GetQ()>0) { | |
741 | d=(AliITSdigitSDD*)digits->UncheckedAt(b->GetIndex()); | |
742 | l0=(d->GetTracks())[0]; | |
743 | } | |
744 | } | |
745 | if (l0<0) { | |
746 | b=&bins[s][idx[k]-nzBins-1]; | |
747 | if (b->GetQ()>0) { | |
748 | d=(AliITSdigitSDD*)digits->UncheckedAt(b->GetIndex()); | |
749 | l0=(d->GetTracks())[0]; | |
750 | } | |
751 | } | |
7f4044f0 | 752 | |
c391f9d9 | 753 | { |
754 | Int_t lab[3]; | |
755 | lab[0]=(d->GetTracks())[0]; | |
756 | lab[1]=(d->GetTracks())[1]; | |
757 | lab[2]=(d->GetTracks())[2]; | |
758 | //CheckLabels(lab); | |
759 | c.SetLabel(lab[0],0); | |
760 | c.SetLabel(lab[1],1); | |
761 | c.SetLabel(lab[2],2); | |
762 | } | |
763 | } | |
7f4044f0 | 764 | |
765 | new (cl[ncl]) AliITSclusterV2(c); ncl++; | |
766 | } | |
767 | } | |
c391f9d9 | 768 | } |
769 | ||
770 | void AliITSclustererV2:: | |
771 | FindClustersSDD(const TClonesArray *digits, TClonesArray *clusters) { | |
772 | //------------------------------------------------------------ | |
773 | // Actual SDD cluster finder | |
774 | //------------------------------------------------------------ | |
775 | Int_t kNzBins = fNzSDD + 2; | |
776 | const Int_t kMAXBIN=kNzBins*(fNySDD+2); | |
777 | ||
778 | AliBin *bins[2]; | |
779 | bins[0]=new AliBin[kMAXBIN]; | |
780 | bins[1]=new AliBin[kMAXBIN]; | |
781 | ||
782 | AliITSdigitSDD *d=0; | |
783 | Int_t i, ndigits=digits->GetEntriesFast(); | |
784 | for (i=0; i<ndigits; i++) { | |
785 | d=(AliITSdigitSDD*)digits->UncheckedAt(i); | |
786 | Int_t y=d->GetCoord2()+1; //y | |
787 | Int_t z=d->GetCoord1()+1; //z | |
788 | Int_t q=d->GetSignal(); | |
789 | if (z <= fNzSDD) { | |
790 | bins[0][y*kNzBins+z].SetQ(q); | |
791 | bins[0][y*kNzBins+z].SetMask(1); | |
792 | bins[0][y*kNzBins+z].SetIndex(i); | |
793 | } else { | |
794 | z-=fNzSDD; | |
795 | bins[1][y*kNzBins+z].SetQ(q); | |
796 | bins[1][y*kNzBins+z].SetMask(1); | |
797 | bins[1][y*kNzBins+z].SetIndex(i); | |
798 | } | |
799 | } | |
800 | ||
801 | FindClustersSDD(bins, kMAXBIN, kNzBins, digits, clusters); | |
7f4044f0 | 802 | |
803 | delete[] bins[0]; | |
804 | delete[] bins[1]; | |
805 | } | |
806 | ||
c391f9d9 | 807 | void AliITSclustererV2::FindClustersSDD(AliITSRawStream* input, |
808 | TClonesArray** clusters) | |
809 | { | |
810 | //------------------------------------------------------------ | |
811 | // Actual SDD cluster finder for raw data | |
812 | //------------------------------------------------------------ | |
813 | Int_t nClustersSDD = 0; | |
814 | Int_t kNzBins = fNzSDD + 2; | |
815 | Int_t kMaxBin = kNzBins * (fNySDD+2); | |
816 | AliBin* bins[2] = {NULL, NULL}; | |
817 | ||
818 | // read raw data input stream | |
819 | while (kTRUE) { | |
820 | Bool_t next = input->Next(); | |
821 | if (!next || input->IsNewModule()) { | |
822 | Int_t iModule = input->GetPrevModuleID(); | |
823 | ||
824 | // when all data from a module was read, search for clusters | |
825 | if (bins[0]) { | |
826 | clusters[iModule] = new TClonesArray("AliITSclusterV2"); | |
827 | fI = iModule; | |
828 | FindClustersSDD(bins, kMaxBin, kNzBins, NULL, clusters[iModule]); | |
829 | Int_t nClusters = clusters[iModule]->GetEntriesFast(); | |
830 | nClustersSDD += nClusters; | |
831 | delete[] bins[0]; | |
832 | delete[] bins[1]; | |
833 | } | |
834 | ||
835 | if (!next) break; | |
836 | bins[0] = new AliBin[kMaxBin]; | |
837 | bins[1] = new AliBin[kMaxBin]; | |
838 | } | |
839 | ||
840 | // fill the current digit into the bins array | |
841 | Int_t iz = input->GetCoord1()+1; | |
842 | Int_t side = ((iz <= fNzSDD) ? 0 : 1); | |
843 | iz -= side*fNzSDD; | |
844 | Int_t index = (input->GetCoord2()+1) * kNzBins + iz; | |
845 | bins[side][index].SetQ(input->GetSignal()); | |
846 | bins[side][index].SetMask(1); | |
847 | bins[side][index].SetIndex(index); | |
848 | } | |
849 | ||
850 | Info("FindClustersSDD", "found clusters in ITS SDD: %d", nClustersSDD); | |
851 | } | |
852 | ||
853 | void AliITSclustererV2:: | |
854 | FindClustersSSD(Ali1Dcluster* neg, Int_t nn, | |
855 | Ali1Dcluster* pos, Int_t np, | |
856 | TClonesArray *clusters) { | |
857 | //------------------------------------------------------------ | |
858 | // Actual SSD cluster finder | |
859 | //------------------------------------------------------------ | |
860 | TClonesArray &cl=*clusters; | |
861 | ||
862 | Int_t lab[4]={-2,-2,-2,-2}; | |
863 | Float_t tanp=fTanP, tann=fTanN; | |
864 | if (fI>fLastSSD1) {tann=fTanP; tanp=fTanN;} | |
865 | ||
866 | Int_t idet=fNdet[fI]; | |
867 | Int_t ncl=0; | |
868 | for (Int_t i=0; i<np; i++) { | |
869 | //Float_t dq_min=1.e+33; | |
870 | Float_t ybest=1000,zbest=1000,qbest=0; | |
871 | Float_t yp=pos[i].GetY()*fYpitchSSD; | |
872 | for (Int_t j=0; j<nn; j++) { | |
873 | //if (pos[i].fTracks[0] != neg[j].fTracks[0]) continue; | |
874 | Float_t yn=neg[j].GetY()*fYpitchSSD; | |
875 | Float_t zt=(2*fHlSSD*tanp + yp - yn)/(tann+tanp); | |
876 | Float_t yt=yn + tann*zt; | |
877 | zt-=fHlSSD; yt-=fHwSSD; | |
878 | if (TMath::Abs(yt)<fHwSSD+0.01) | |
879 | if (TMath::Abs(zt)<fHlSSD+0.01) { | |
880 | //if (TMath::Abs(pos[i].GetQ()-neg[j].GetQ())<dq_min) { | |
881 | //dq_min=TMath::Abs(pos[i].GetQ()-neg[j].GetQ()); | |
882 | ybest=yt; zbest=zt; | |
883 | qbest=0.5*(pos[i].GetQ()+neg[j].GetQ()); | |
884 | ||
885 | lab[0]=pos[i].GetLabel(0); | |
886 | lab[1]=pos[i].GetLabel(1); | |
887 | lab[2]=neg[i].GetLabel(0); | |
888 | lab[3]=(((i<<10) + j)<<10) + idet; // pos|neg|det | |
889 | Float_t lp[5]; | |
890 | lp[0]=-ybest-fYshift[fI]; | |
891 | lp[1]= zbest+fZshift[fI]; | |
892 | lp[2]=0.0025*0.0025; //SigmaY2 | |
893 | lp[3]=0.110*0.110; //SigmaZ2 | |
894 | if (pos[i].GetNd()+neg[j].GetNd() > 4) { | |
895 | lp[2]*=9; | |
896 | lp[3]*=9; | |
897 | } | |
898 | lp[4]=qbest; //Q | |
899 | ||
900 | //CheckLabels(lab); | |
901 | new (cl[ncl]) AliITSclusterV2(lab,lp); ncl++; | |
902 | } | |
903 | } | |
904 | /* | |
905 | if (ybest<100) { | |
906 | lab[3]=idet; | |
907 | Float_t lp[5]; | |
908 | lp[0]=-ybest-fYshift[fI]; | |
909 | lp[1]= zbest+fZshift[fI]; | |
910 | lp[2]=0.002*0.002; //SigmaY2 | |
911 | lp[3]=0.080*0.080; //SigmaZ2 | |
912 | lp[4]=qbest; //Q | |
913 | // | |
914 | new (cl[ncl]) AliITSclusterV2(lab,lp); ncl++; | |
915 | } | |
916 | */ | |
917 | } | |
918 | } | |
919 | ||
7f4044f0 | 920 | void AliITSclustererV2:: |
921 | FindClustersSSD(const TClonesArray *digits, TClonesArray *clusters) { | |
922 | //------------------------------------------------------------ | |
923 | // Actual SSD cluster finder | |
924 | //------------------------------------------------------------ | |
925 | Int_t smax=digits->GetEntriesFast(); | |
926 | if (smax==0) return; | |
927 | ||
928 | const Int_t MAX=1000; | |
929 | Int_t np=0, nn=0; | |
930 | Ali1Dcluster pos[MAX], neg[MAX]; | |
931 | Float_t y=0., q=0., qmax=0.; | |
932 | Int_t lab[4]={-2,-2,-2,-2}; | |
933 | ||
7f4044f0 | 934 | AliITSdigitSSD *d=(AliITSdigitSSD*)digits->UncheckedAt(0); |
935 | q += d->GetSignal(); | |
936 | y += d->GetCoord2()*d->GetSignal(); | |
937 | qmax=d->GetSignal(); | |
938 | lab[0]=d->GetTrack(0); lab[1]=d->GetTrack(1); lab[2]=d->GetTrack(2); | |
939 | Int_t curr=d->GetCoord2(); | |
940 | Int_t flag=d->GetCoord1(); | |
941 | Int_t *n=&nn; | |
942 | Ali1Dcluster *c=neg; | |
c391f9d9 | 943 | Int_t nd=1; |
7f4044f0 | 944 | for (Int_t s=1; s<smax; s++) { |
945 | d=(AliITSdigitSSD*)digits->UncheckedAt(s); | |
946 | Int_t strip=d->GetCoord2(); | |
947 | if ((strip-curr) > 1 || flag!=d->GetCoord1()) { | |
948 | c[*n].SetY(y/q); | |
949 | c[*n].SetQ(q); | |
950 | c[*n].SetNd(nd); | |
951 | c[*n].SetLabels(lab); | |
952 | //Split suspiciously big cluster | |
953 | if (nd>3) { | |
954 | c[*n].SetY(y/q-0.5*nd); | |
955 | c[*n].SetQ(0.5*q); | |
956 | (*n)++; | |
957 | if (*n==MAX) { | |
958 | cerr<< | |
959 | "AliITSclustererV2::FindClustersSSD: Too many 1D clusters !\n"; | |
960 | return; | |
961 | } | |
962 | c[*n].SetY(y/q+0.5*nd); | |
963 | c[*n].SetQ(0.5*q); | |
964 | c[*n].SetNd(nd); | |
965 | c[*n].SetLabels(lab); | |
966 | } | |
967 | (*n)++; | |
968 | if (*n==MAX) { | |
969 | cerr<<"AliITSclustererV2::FindClustersSSD: Too many 1D clusters !\n"; | |
970 | return; | |
971 | } | |
972 | y=q=qmax=0.; | |
973 | nd=0; | |
974 | lab[0]=lab[1]=lab[2]=-2; | |
975 | if (flag!=d->GetCoord1()) { n=&np; c=pos; } | |
976 | } | |
977 | flag=d->GetCoord1(); | |
978 | q += d->GetSignal(); | |
979 | y += d->GetCoord2()*d->GetSignal(); | |
980 | nd++; | |
981 | if (d->GetSignal()>qmax) { | |
982 | qmax=d->GetSignal(); | |
983 | lab[0]=d->GetTrack(0); lab[1]=d->GetTrack(1); lab[2]=d->GetTrack(2); | |
984 | } | |
985 | curr=strip; | |
986 | } | |
987 | c[*n].SetY(y/q); | |
988 | c[*n].SetQ(q); | |
989 | c[*n].SetNd(nd); | |
990 | c[*n].SetLabels(lab); | |
991 | //Split suspiciously big cluster | |
992 | if (nd>3) { | |
993 | c[*n].SetY(y/q-0.5*nd); | |
994 | c[*n].SetQ(0.5*q); | |
995 | (*n)++; | |
996 | if (*n==MAX) { | |
997 | cerr<<"AliITSclustererV2::FindClustersSSD: Too many 1D clusters !\n"; | |
998 | return; | |
999 | } | |
1000 | c[*n].SetY(y/q+0.5*nd); | |
1001 | c[*n].SetQ(0.5*q); | |
1002 | c[*n].SetNd(nd); | |
1003 | c[*n].SetLabels(lab); | |
1004 | } | |
1005 | (*n)++; | |
1006 | if (*n==MAX) { | |
1007 | cerr<<"AliITSclustererV2::FindClustersSSD: Too many 1D clusters !\n"; | |
1008 | return; | |
1009 | } | |
1010 | ||
c391f9d9 | 1011 | FindClustersSSD(neg, nn, pos, np, clusters); |
1012 | } | |
7f4044f0 | 1013 | |
c391f9d9 | 1014 | void AliITSclustererV2::FindClustersSSD(AliITSRawStream* input, |
1015 | TClonesArray** clusters) | |
1016 | { | |
1017 | //------------------------------------------------------------ | |
1018 | // Actual SSD cluster finder for raw data | |
1019 | //------------------------------------------------------------ | |
1020 | Int_t nClustersSSD = 0; | |
1021 | const Int_t MAX = 1000; | |
1022 | Ali1Dcluster clusters1D[2][MAX]; | |
1023 | Int_t nClusters[2] = {0, 0}; | |
1024 | Float_t q = 0.; | |
1025 | Float_t y = 0.; | |
1026 | Int_t nDigits = 0; | |
1027 | Int_t prevStrip = -1; | |
1028 | Int_t prevFlag = -1; | |
1029 | ||
1030 | // read raw data input stream | |
1031 | while (kTRUE) { | |
1032 | Bool_t next = input->Next(); | |
1033 | ||
1034 | // check if a new cluster starts | |
1035 | Int_t strip = input->GetCoord2(); | |
1036 | Int_t flag = input->GetCoord1(); | |
1037 | if ((!next || input->IsNewModule() || | |
1038 | (strip-prevStrip > 1) || (flag != prevFlag)) && | |
1039 | (nDigits > 0)) { | |
1040 | if (nClusters[prevFlag] == MAX) { | |
1041 | Error("FindClustersSSD", "Too many 1D clusters !"); | |
1042 | return; | |
1043 | } | |
1044 | Ali1Dcluster& cluster = clusters1D[prevFlag][nClusters[prevFlag]++]; | |
1045 | cluster.SetY(y/q); | |
1046 | cluster.SetQ(q); | |
1047 | cluster.SetNd(nDigits); | |
1048 | ||
1049 | //Split suspiciously big cluster | |
1050 | if (nDigits > 3) { | |
1051 | cluster.SetY(y/q - 0.5*nDigits); | |
1052 | cluster.SetQ(0.5*q); | |
1053 | if (nClusters[prevFlag] == MAX) { | |
1054 | Error("FindClustersSSD", "Too many 1D clusters !"); | |
1055 | return; | |
1056 | } | |
1057 | Ali1Dcluster& cluster2 = clusters1D[prevFlag][nClusters[prevFlag]++]; | |
1058 | cluster2.SetY(y/q + 0.5*nDigits); | |
1059 | cluster2.SetQ(0.5*q); | |
1060 | cluster2.SetNd(nDigits); | |
7f4044f0 | 1061 | } |
c391f9d9 | 1062 | y = q = 0.; |
1063 | nDigits = 0; | |
7f4044f0 | 1064 | } |
c391f9d9 | 1065 | |
1066 | if (!next || input->IsNewModule()) { | |
1067 | Int_t iModule = input->GetPrevModuleID(); | |
1068 | ||
1069 | // when all data from a module was read, search for clusters | |
1070 | if (prevFlag >= 0) { | |
1071 | clusters[iModule] = new TClonesArray("AliITSclusterV2"); | |
1072 | fI = iModule; | |
1073 | FindClustersSSD(&clusters1D[0][0], nClusters[0], | |
1074 | &clusters1D[1][0], nClusters[1], clusters[iModule]); | |
1075 | Int_t nClusters = clusters[iModule]->GetEntriesFast(); | |
1076 | nClustersSSD += nClusters; | |
1077 | } | |
1078 | ||
1079 | if (!next) break; | |
1080 | nClusters[0] = nClusters[1] = 0; | |
1081 | y = q = 0.; | |
1082 | nDigits = 0; | |
7f4044f0 | 1083 | } |
c391f9d9 | 1084 | |
1085 | // add digit to current cluster | |
1086 | q += input->GetSignal(); | |
1087 | y += strip * input->GetSignal(); | |
1088 | nDigits++; | |
1089 | prevStrip = strip; | |
1090 | prevFlag = flag; | |
7f4044f0 | 1091 | } |
c391f9d9 | 1092 | |
1093 | Info("FindClustersSSD", "found clusters in ITS SSD: %d", nClustersSSD); | |
7f4044f0 | 1094 | } |
1095 | ||
1096 | #else //V1 | |
1097 | ||
1098 | #include "AliITSDetType.h" | |
1099 | ||
1100 | #include "AliITSsegmentationSPD.h" | |
1101 | #include "AliITSClusterFinderSPD.h" | |
1102 | ||
1103 | #include "AliITSresponseSDD.h" | |
1104 | #include "AliITSsegmentationSDD.h" | |
1105 | #include "AliITSClusterFinderSDD.h" | |
1106 | ||
1107 | #include "AliITSsegmentationSSD.h" | |
1108 | #include "AliITSClusterFinderSSD.h" | |
1109 | ||
1110 | ||
1111 | void AliITSclustererV2:: | |
1112 | FindClustersSPD(const TClonesArray *digits, TClonesArray *clusters) { | |
1113 | //------------------------------------------------------------ | |
1114 | // Actual SPD cluster finding based on AliITSClusterFinderSPD | |
1115 | //------------------------------------------------------------ | |
1116 | static AliITS *its=(AliITS*)gAlice->GetModule("ITS"); | |
1117 | static TClonesArray *points=its->RecPoints(); | |
1118 | static AliITSsegmentationSPD *seg= | |
1119 | (AliITSsegmentationSPD *)its->DetType(0)->GetSegmentationModel(); | |
1120 | static AliITSClusterFinderSPD cf(seg, (TClonesArray*)digits, points); | |
1121 | ||
1122 | cf.FindRawClusters(fI); | |
1cca57bf | 1123 | RecPoints2Clusters(points, fI, clusters); |
7f4044f0 | 1124 | its->ResetRecPoints(); |
1125 | ||
1126 | } | |
1127 | ||
1128 | void AliITSclustererV2:: | |
1129 | FindClustersSDD(const TClonesArray *digits, TClonesArray *clusters) { | |
1130 | //------------------------------------------------------------ | |
1131 | // Actual SDD cluster finding based on AliITSClusterFinderSDD | |
1132 | //------------------------------------------------------------ | |
1133 | static AliITS *its=(AliITS*)gAlice->GetModule("ITS"); | |
1134 | static TClonesArray *points=its->RecPoints(); | |
1135 | static AliITSresponseSDD *resp= | |
1136 | (AliITSresponseSDD *)its->DetType(1)->GetResponseModel(); | |
1137 | static AliITSsegmentationSDD *seg= | |
1138 | (AliITSsegmentationSDD *)its->DetType(1)->GetSegmentationModel(); | |
1139 | static AliITSClusterFinderSDD | |
1140 | cf(seg,resp,(TClonesArray*)digits,its->ClustersAddress(1)); | |
1141 | ||
1142 | cf.FindRawClusters(fI); | |
1cca57bf | 1143 | Int_t nc=points->GetEntriesFast(); |
1144 | while (nc--) { //To be consistent with the SSD cluster charges | |
1145 | AliITSRecPoint *p=(AliITSRecPoint*)points->UncheckedAt(nc); | |
1146 | p->SetQ(p->GetQ/12.); | |
1147 | } | |
1148 | RecPoints2Clusters(points, fI, clusters); | |
7f4044f0 | 1149 | its->ResetClusters(1); |
1150 | its->ResetRecPoints(); | |
1151 | ||
1152 | } | |
1153 | ||
1154 | void AliITSclustererV2:: | |
1155 | FindClustersSSD(const TClonesArray *digits, TClonesArray *clusters) { | |
1156 | //------------------------------------------------------------ | |
1157 | // Actual SSD cluster finding based on AliITSClusterFinderSSD | |
1158 | //------------------------------------------------------------ | |
1159 | static AliITS *its=(AliITS*)gAlice->GetModule("ITS"); | |
1160 | static TClonesArray *points=its->RecPoints(); | |
1161 | static AliITSsegmentationSSD *seg= | |
1162 | (AliITSsegmentationSSD *)its->DetType(2)->GetSegmentationModel(); | |
1163 | static AliITSClusterFinderSSD cf(seg,(TClonesArray*)digits); | |
1164 | ||
1165 | cf.FindRawClusters(fI); | |
1cca57bf | 1166 | RecPoints2Clusters(points, fI, clusters); |
7f4044f0 | 1167 | its->ResetRecPoints(); |
1168 | ||
1169 | } | |
1170 | ||
1171 | #endif |