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