]>
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); | |
e43c066c | 54 | fNlayer[m] = lay-1; |
7f4044f0 | 55 | } |
c391f9d9 | 56 | fNModules = g->GetIndexMax(); |
7f4044f0 | 57 | |
58 | //SPD geometry | |
59 | fLastSPD1=g->GetModuleIndex(2,1,1)-1; | |
60 | fNySPD=256; fNzSPD=160; | |
61 | fYpitchSPD=0.0050; | |
62 | fZ1pitchSPD=0.0425; fZ2pitchSPD=0.0625; | |
63 | fHwSPD=0.64; fHlSPD=3.48; | |
64 | fYSPD[0]=0.5*fYpitchSPD; | |
65 | for (m=1; m<fNySPD; m++) fYSPD[m]=fYSPD[m-1]+fYpitchSPD; | |
66 | fZSPD[0]=fZ1pitchSPD; | |
67 | for (m=1; m<fNzSPD; m++) { | |
68 | Double_t dz=fZ1pitchSPD; | |
69 | if (m==31 || m==32 || m==63 || m==64 || m==95 || m==96 || | |
70 | m==127 || m==128) dz=fZ2pitchSPD; | |
71 | fZSPD[m]=fZSPD[m-1]+dz; | |
72 | } | |
73 | for (m=0; m<fNzSPD; m++) { | |
74 | Double_t dz=0.5*fZ1pitchSPD; | |
75 | if (m==31 || m==32 || m==63 || m==64 || m==95 || m==96 || | |
76 | m==127 || m==128) dz=0.5*fZ2pitchSPD; | |
77 | fZSPD[m]-=dz; | |
78 | } | |
79 | ||
80 | //SDD geometry | |
81 | fNySDD=256; fNzSDD=256; | |
82 | fYpitchSDD=0.01825; | |
83 | fZpitchSDD=0.02940; | |
84 | fHwSDD=3.5085; fHlSDD=3.7632; | |
85 | fYoffSDD=0.0425; | |
86 | ||
87 | //SSD geometry | |
88 | fLastSSD1=g->GetModuleIndex(6,1,1)-1; | |
89 | fYpitchSSD=0.0095; | |
90 | fHwSSD=3.65; | |
91 | fHlSSD=2.00; | |
92 | fTanP=0.0275; | |
93 | fTanN=0.0075; | |
5102d526 | 94 | |
7f4044f0 | 95 | } |
96 | ||
e3d91d99 | 97 | |
c630aafd | 98 | Int_t AliITSclustererV2::Digits2Clusters(TTree *dTree, TTree *cTree) { |
7f4044f0 | 99 | //------------------------------------------------------------ |
100 | // This function creates ITS clusters | |
101 | //------------------------------------------------------------ | |
102 | Int_t ncl=0; | |
7f4044f0 | 103 | |
7f4044f0 | 104 | if (!dTree) { |
c630aafd | 105 | Error("Digits2Clusters","Can't get the tree with digits !"); |
106 | return 1; | |
7f4044f0 | 107 | } |
108 | ||
109 | TClonesArray *digitsSPD=new TClonesArray("AliITSdigitSPD",3000); | |
110 | dTree->SetBranchAddress("ITSDigitsSPD",&digitsSPD); | |
111 | TClonesArray *digitsSDD=new TClonesArray("AliITSdigitSDD",3000); | |
112 | dTree->SetBranchAddress("ITSDigitsSDD",&digitsSDD); | |
113 | TClonesArray *digitsSSD=new TClonesArray("AliITSdigitSSD",3000); | |
114 | dTree->SetBranchAddress("ITSDigitsSSD",&digitsSSD); | |
115 | ||
7f4044f0 | 116 | TClonesArray *clusters=new TClonesArray("AliITSclusterV2",1000); |
c630aafd | 117 | TBranch *branch=cTree->GetBranch("Clusters"); |
118 | if (!branch) cTree->Branch("Clusters",&clusters); | |
119 | else branch->SetAddress(&clusters); | |
120 | ||
121 | Int_t mmax=(Int_t)dTree->GetEntries(); | |
cbca22ab | 122 | if (mmax!=fNModules) { |
123 | Error("Digits2Clusters","Number of entries != number of modules !"); | |
124 | return 1; | |
125 | } | |
7f4044f0 | 126 | |
127 | for (fI=0; fI<mmax; fI++) { | |
128 | dTree->GetEvent(fI); | |
129 | ||
130 | if (digitsSPD->GetEntriesFast()!=0) | |
babd135a | 131 | FindClustersSPD(digitsSPD,clusters); |
132 | else | |
133 | if(digitsSDD->GetEntriesFast()!=0) | |
134 | FindClustersSDD(digitsSDD,clusters); | |
135 | else if(digitsSSD->GetEntriesFast()!=0) | |
136 | FindClustersSSD(digitsSSD,clusters); | |
137 | ||
7f4044f0 | 138 | ncl+=clusters->GetEntriesFast(); |
139 | ||
c630aafd | 140 | cTree->Fill(); |
7f4044f0 | 141 | |
142 | digitsSPD->Clear(); | |
143 | digitsSDD->Clear(); | |
144 | digitsSSD->Clear(); | |
145 | clusters->Clear(); | |
146 | } | |
c630aafd | 147 | |
148 | //cTree->Write(); | |
7f4044f0 | 149 | |
150 | delete clusters; | |
151 | ||
152 | delete digitsSPD; | |
153 | delete digitsSDD; | |
154 | delete digitsSSD; | |
155 | ||
7941072e | 156 | //delete dTree; |
157 | ||
c630aafd | 158 | Info("Digits2Clusters","Number of found clusters : %d",ncl); |
7f4044f0 | 159 | |
c630aafd | 160 | return 0; |
7f4044f0 | 161 | } |
162 | ||
7941072e | 163 | void AliITSclustererV2::Digits2Clusters(AliRawReader* rawReader) { |
c391f9d9 | 164 | //------------------------------------------------------------ |
165 | // This function creates ITS clusters from raw data | |
166 | //------------------------------------------------------------ | |
7941072e | 167 | AliRunLoader* runLoader = AliRunLoader::GetRunLoader(); |
168 | if (!runLoader) { | |
169 | Error("Digits2Clusters", "no run loader found"); | |
170 | return; | |
171 | } | |
babd135a | 172 | runLoader->LoadKinematics(); |
7941072e | 173 | AliLoader* itsLoader = runLoader->GetLoader("ITSLoader"); |
174 | if (!itsLoader) { | |
175 | Error("Digits2Clusters", "no loader for ITS found"); | |
c391f9d9 | 176 | return; |
177 | } | |
7941072e | 178 | if (!itsLoader->TreeR()) itsLoader->MakeTree("R"); |
179 | TTree* cTree = itsLoader->TreeR(); | |
c391f9d9 | 180 | |
c391f9d9 | 181 | TClonesArray *array=new TClonesArray("AliITSclusterV2",1000); |
7941072e | 182 | cTree->Branch("Clusters",&array); |
c391f9d9 | 183 | delete array; |
184 | ||
f1e0c1c5 | 185 | TClonesArray** clusters = new TClonesArray*[fNModules]; |
04fa961a | 186 | for (Int_t iModule = 0; iModule < fNModules; iModule++) { |
187 | clusters[iModule] = NULL; | |
188 | } | |
c391f9d9 | 189 | // one TClonesArray per module |
190 | ||
7941072e | 191 | rawReader->Reset(); |
192 | AliITSRawStreamSPD inputSPD(rawReader); | |
c391f9d9 | 193 | FindClustersSPD(&inputSPD, clusters); |
7941072e | 194 | |
195 | rawReader->Reset(); | |
196 | AliITSRawStreamSDD inputSDD(rawReader); | |
c391f9d9 | 197 | FindClustersSDD(&inputSDD, clusters); |
7941072e | 198 | |
199 | rawReader->Reset(); | |
200 | AliITSRawStreamSSD inputSSD(rawReader); | |
c391f9d9 | 201 | FindClustersSSD(&inputSSD, clusters); |
202 | ||
203 | // write all clusters to the tree | |
204 | Int_t nClusters = 0; | |
205 | for (Int_t iModule = 0; iModule < fNModules; iModule++) { | |
04fa961a | 206 | array = clusters[iModule]; |
c391f9d9 | 207 | if (!array) { |
208 | Error("Digits2Clusters", "data for module %d missing!", iModule); | |
209 | array = new TClonesArray("AliITSclusterV2"); | |
210 | } | |
7941072e | 211 | cTree->SetBranchAddress("Clusters", &array); |
212 | cTree->Fill(); | |
c391f9d9 | 213 | nClusters += array->GetEntriesFast(); |
214 | delete array; | |
215 | } | |
7941072e | 216 | itsLoader->WriteRecPoints("OVERWRITE"); |
c391f9d9 | 217 | |
04fa961a | 218 | delete[] clusters; |
219 | ||
c391f9d9 | 220 | Info("Digits2Clusters", "total number of found clusters in ITS: %d\n", |
221 | nClusters); | |
222 | } | |
223 | ||
7f4044f0 | 224 | //**** Fast clusters ******************************* |
225 | #include "TParticle.h" | |
226 | ||
227 | #include "AliITS.h" | |
228 | #include "AliITSmodule.h" | |
229 | #include "AliITSRecPoint.h" | |
230 | #include "AliITSsimulationFastPoints.h" | |
231 | #include "AliITSRecPoint.h" | |
232 | ||
babd135a | 233 | /* |
7f4044f0 | 234 | static void CheckLabels(Int_t lab[3]) { |
235 | //------------------------------------------------------------ | |
236 | // Tries to find mother's labels | |
237 | //------------------------------------------------------------ | |
238 | Int_t label=lab[0]; | |
239 | if (label>=0) { | |
5d12ce38 | 240 | TParticle *part=(TParticle*)gAlice->GetMCApp()->Particle(label); |
7f4044f0 | 241 | label=-3; |
242 | while (part->P() < 0.005) { | |
243 | Int_t m=part->GetFirstMother(); | |
c630aafd | 244 | if (m<0) { |
245 | Info("CheckLabels","Primary momentum: %f",part->P()); | |
246 | break; | |
247 | } | |
1cca57bf | 248 | if (part->GetStatusCode()>0) { |
c630aafd | 249 | Info("CheckLabels","Primary momentum: %f",part->P()); |
250 | break; | |
1cca57bf | 251 | } |
7f4044f0 | 252 | label=m; |
5d12ce38 | 253 | part=(TParticle*)gAlice->GetMCApp()->Particle(label); |
7f4044f0 | 254 | } |
255 | if (lab[1]<0) lab[1]=label; | |
256 | else if (lab[2]<0) lab[2]=label; | |
257 | else ;//cerr<<"CheckLabels : No empty labels !\n"; | |
258 | } | |
259 | } | |
babd135a | 260 | */ |
261 | static void CheckLabels(Int_t lab[3]) { | |
262 | //------------------------------------------------------------ | |
263 | // Tries to find mother's labels | |
264 | //------------------------------------------------------------ | |
265 | ||
0e986ef9 | 266 | if(lab[0]<0 && lab[1]<0 && lab[2]<0) return; // In case of no labels just exit |
267 | ||
babd135a | 268 | Int_t ntracks = gAlice->GetMCApp()->GetNtrack(); |
269 | for (Int_t i=0;i<3;i++){ | |
270 | Int_t label = lab[i]; | |
271 | if (label>=0 && label<ntracks) { | |
272 | TParticle *part=(TParticle*)gAlice->GetMCApp()->Particle(label); | |
273 | if (part->P() < 0.005) { | |
274 | Int_t m=part->GetFirstMother(); | |
275 | if (m<0) { | |
276 | continue; | |
277 | } | |
278 | if (part->GetStatusCode()>0) { | |
279 | continue; | |
280 | } | |
281 | lab[i]=m; | |
282 | } | |
283 | } | |
284 | } | |
285 | ||
286 | } | |
287 | ||
288 | static void CheckLabels2(Int_t lab[10]) { | |
289 | //------------------------------------------------------------ | |
290 | // Tries to find mother's labels | |
291 | //------------------------------------------------------------ | |
e43c066c | 292 | Int_t nlabels =0; |
293 | for (Int_t i=0;i<10;i++) if (lab[i]>=0) nlabels++; | |
0e986ef9 | 294 | if(nlabels == 0) return; // In case of no labels just exit |
295 | ||
296 | Int_t ntracks = gAlice->GetMCApp()->GetNtrack(); | |
e43c066c | 297 | |
babd135a | 298 | for (Int_t i=0;i<10;i++){ |
299 | Int_t label = lab[i]; | |
300 | if (label>=0 && label<ntracks) { | |
301 | TParticle *part=(TParticle*)gAlice->GetMCApp()->Particle(label); | |
e43c066c | 302 | if (part->P() < 0.02) { |
303 | Int_t m=part->GetFirstMother(); | |
304 | if (m<0) { | |
305 | continue; | |
306 | } | |
307 | if (part->GetStatusCode()>0) { | |
308 | continue; | |
309 | } | |
310 | lab[i]=m; | |
311 | } | |
312 | else | |
313 | if (part->P() < 0.12 && nlabels>3) { | |
314 | lab[i]=-2; | |
315 | nlabels--; | |
316 | } | |
317 | } | |
318 | else{ | |
319 | if ( (label>ntracks||label <0) && nlabels>3) { | |
320 | lab[i]=-2; | |
321 | nlabels--; | |
322 | } | |
323 | } | |
324 | } | |
325 | if (nlabels>3){ | |
326 | for (Int_t i=0;i<10;i++){ | |
327 | if (nlabels>3){ | |
328 | Int_t label = lab[i]; | |
329 | if (label>=0 && label<ntracks) { | |
330 | TParticle *part=(TParticle*)gAlice->GetMCApp()->Particle(label); | |
331 | if (part->P() < 0.1) { | |
332 | lab[i]=-2; | |
333 | nlabels--; | |
334 | } | |
babd135a | 335 | } |
babd135a | 336 | } |
e43c066c | 337 | } |
babd135a | 338 | } |
e43c066c | 339 | |
babd135a | 340 | //compress labels -- if multi-times the same |
341 | Int_t lab2[10]; | |
342 | for (Int_t i=0;i<10;i++) lab2[i]=-2; | |
343 | for (Int_t i=0;i<10 ;i++){ | |
344 | if (lab[i]<0) continue; | |
345 | for (Int_t j=0;j<10 &&lab2[j]!=lab[i];j++){ | |
346 | if (lab2[j]<0) { | |
347 | lab2[j]= lab[i]; | |
348 | break; | |
349 | } | |
350 | } | |
351 | } | |
352 | for (Int_t j=0;j<10;j++) lab[j]=lab2[j]; | |
353 | ||
354 | } | |
355 | ||
356 | static void AddLabel(Int_t lab[10], Int_t label) { | |
0e986ef9 | 357 | |
358 | if(label<0) return; // In case of no label just exit | |
359 | ||
e43c066c | 360 | Int_t ntracks = gAlice->GetMCApp()->GetNtrack(); |
361 | if (label>ntracks) return; | |
babd135a | 362 | for (Int_t i=0;i<10;i++){ |
0e986ef9 | 363 | // if (label<0) break; |
babd135a | 364 | if (lab[i]==label) break; |
365 | if (lab[i]<0) { | |
366 | lab[i]= label; | |
367 | break; | |
368 | } | |
369 | } | |
370 | } | |
7f4044f0 | 371 | |
1cca57bf | 372 | void AliITSclustererV2::RecPoints2Clusters |
373 | (const TClonesArray *points, Int_t idx, TClonesArray *clusters) { | |
7f4044f0 | 374 | //------------------------------------------------------------ |
1cca57bf | 375 | // Conversion AliITSRecPoint -> AliITSclusterV2 for the ITS |
376 | // subdetector indexed by idx | |
7f4044f0 | 377 | //------------------------------------------------------------ |
378 | TClonesArray &cl=*clusters; | |
379 | Int_t ncl=points->GetEntriesFast(); | |
380 | for (Int_t i=0; i<ncl; i++) { | |
381 | AliITSRecPoint *p = (AliITSRecPoint *)points->UncheckedAt(i); | |
382 | Float_t lp[5]; | |
f363d377 | 383 | lp[0]=-(-p->GetX()+fYshift[idx]); if (idx<=fLastSPD1) lp[0]*=-1; //SPD1 |
384 | lp[1]= -p->GetZ()+fZshift[idx]; | |
1cca57bf | 385 | lp[2]=p->GetSigmaX2(); |
386 | lp[3]=p->GetSigmaZ2(); | |
b6087704 | 387 | lp[4]=p->GetQ()*36./23333.; //electrons -> ADC |
7f4044f0 | 388 | Int_t lab[4]; |
389 | lab[0]=p->GetLabel(0); lab[1]=p->GetLabel(1); lab[2]=p->GetLabel(2); | |
1cca57bf | 390 | lab[3]=fNdet[idx]; |
7f4044f0 | 391 | CheckLabels(lab); |
babd135a | 392 | Int_t dummy[3]={0,0,0}; |
393 | new (cl[i]) AliITSclusterV2(lab,lp, dummy); | |
7f4044f0 | 394 | } |
395 | } | |
396 | ||
c630aafd | 397 | Int_t AliITSclustererV2::Hits2Clusters(TTree *hTree, TTree *cTree) { |
7f4044f0 | 398 | //------------------------------------------------------------ |
399 | // This function creates ITS clusters | |
400 | //------------------------------------------------------------ | |
7f4044f0 | 401 | if (!gAlice) { |
c630aafd | 402 | Error("Hits2Clusters","gAlice==0 !"); |
403 | return 1; | |
7f4044f0 | 404 | } |
405 | ||
406 | AliITS *its = (AliITS*)gAlice->GetModule("ITS"); | |
407 | if (!its) { | |
c630aafd | 408 | Error("Hits2Clusters","Can't find the ITS !"); |
409 | return 2; | |
7f4044f0 | 410 | } |
411 | AliITSgeom *geom=its->GetITSgeom(); | |
412 | Int_t mmax=geom->GetIndexMax(); | |
413 | ||
414 | its->InitModules(-1,mmax); | |
c630aafd | 415 | its->FillModules(hTree,0); |
7f4044f0 | 416 | |
7f4044f0 | 417 | TClonesArray *clusters=new TClonesArray("AliITSclusterV2",1000); |
c630aafd | 418 | TBranch *branch=cTree->GetBranch("Clusters"); |
419 | if (!branch) cTree->Branch("Clusters",&clusters); | |
420 | else branch->SetAddress(&clusters); | |
7f4044f0 | 421 | |
422 | static TClonesArray *points=its->RecPoints(); | |
423 | AliITSsimulationFastPoints sim; | |
424 | Int_t ncl=0; | |
425 | for (Int_t m=0; m<mmax; m++) { | |
426 | AliITSmodule *mod=its->GetModule(m); | |
427 | sim.CreateFastRecPoints(mod,m,gRandom); | |
428 | ||
1cca57bf | 429 | RecPoints2Clusters(points, m, clusters); |
7f4044f0 | 430 | its->ResetRecPoints(); |
431 | ||
432 | ncl+=clusters->GetEntriesFast(); | |
c630aafd | 433 | cTree->Fill(); |
7f4044f0 | 434 | clusters->Clear(); |
435 | } | |
7f4044f0 | 436 | |
c630aafd | 437 | Info("Hits2Clusters","Number of found fast clusters : %d",ncl); |
438 | ||
439 | //cTree->Write(); | |
7f4044f0 | 440 | |
441 | delete clusters; | |
442 | ||
c630aafd | 443 | return 0; |
7f4044f0 | 444 | } |
445 | ||
446 | //*********************************** | |
447 | ||
448 | #ifndef V1 | |
449 | ||
450 | void AliITSclustererV2:: | |
451 | FindCluster(Int_t k,Int_t maxz,AliBin *bins,Int_t &n,Int_t *idx) { | |
452 | //------------------------------------------------------------ | |
453 | // returns an array of indices of digits belonging to the cluster | |
454 | // (needed when the segmentation is not regular) | |
455 | //------------------------------------------------------------ | |
456 | if (n<200) idx[n++]=bins[k].GetIndex(); | |
457 | bins[k].Use(); | |
458 | ||
459 | if (bins[k-maxz].IsNotUsed()) FindCluster(k-maxz,maxz,bins,n,idx); | |
460 | if (bins[k-1 ].IsNotUsed()) FindCluster(k-1 ,maxz,bins,n,idx); | |
461 | if (bins[k+maxz].IsNotUsed()) FindCluster(k+maxz,maxz,bins,n,idx); | |
462 | if (bins[k+1 ].IsNotUsed()) FindCluster(k+1 ,maxz,bins,n,idx); | |
463 | /* | |
464 | if (bins[k-maxz-1].IsNotUsed()) FindCluster(k-maxz-1,maxz,bins,n,idx); | |
465 | if (bins[k-maxz+1].IsNotUsed()) FindCluster(k-maxz+1,maxz,bins,n,idx); | |
466 | if (bins[k+maxz-1].IsNotUsed()) FindCluster(k+maxz-1,maxz,bins,n,idx); | |
467 | if (bins[k+maxz+1].IsNotUsed()) FindCluster(k+maxz+1,maxz,bins,n,idx); | |
468 | */ | |
469 | } | |
470 | ||
471 | void AliITSclustererV2:: | |
472 | FindClustersSPD(const TClonesArray *digits, TClonesArray *clusters) { | |
473 | //------------------------------------------------------------ | |
474 | // Actual SPD cluster finder | |
475 | //------------------------------------------------------------ | |
c391f9d9 | 476 | Int_t kNzBins = fNzSPD + 2; |
477 | const Int_t kMAXBIN=kNzBins*(fNySPD+2); | |
7f4044f0 | 478 | |
479 | Int_t ndigits=digits->GetEntriesFast(); | |
480 | AliBin *bins=new AliBin[kMAXBIN]; | |
481 | ||
482 | Int_t k; | |
483 | AliITSdigitSPD *d=0; | |
484 | for (k=0; k<ndigits; k++) { | |
485 | d=(AliITSdigitSPD*)digits->UncheckedAt(k); | |
486 | Int_t i=d->GetCoord2()+1; //y | |
487 | Int_t j=d->GetCoord1()+1; | |
c391f9d9 | 488 | bins[i*kNzBins+j].SetIndex(k); |
489 | bins[i*kNzBins+j].SetMask(1); | |
7f4044f0 | 490 | } |
491 | ||
492 | Int_t n=0; TClonesArray &cl=*clusters; | |
493 | for (k=0; k<kMAXBIN; k++) { | |
494 | if (!bins[k].IsNotUsed()) continue; | |
495 | Int_t ni=0, idx[200]; | |
c391f9d9 | 496 | FindCluster(k,kNzBins,bins,ni,idx); |
c630aafd | 497 | if (ni==200) { |
498 | Info("FindClustersSPD","Too big cluster !"); | |
499 | continue; | |
500 | } | |
babd135a | 501 | Int_t milab[10]; |
502 | for (Int_t ilab=0;ilab<10;ilab++){ | |
503 | milab[ilab]=-2; | |
504 | } | |
505 | ||
7f4044f0 | 506 | d=(AliITSdigitSPD*)digits->UncheckedAt(idx[0]); |
507 | Int_t ymin=d->GetCoord2(),ymax=ymin; | |
508 | Int_t zmin=d->GetCoord1(),zmax=zmin; | |
e43c066c | 509 | |
7f4044f0 | 510 | for (Int_t l=0; l<ni; l++) { |
511 | d=(AliITSdigitSPD*)digits->UncheckedAt(idx[l]); | |
512 | ||
513 | if (ymin > d->GetCoord2()) ymin=d->GetCoord2(); | |
514 | if (ymax < d->GetCoord2()) ymax=d->GetCoord2(); | |
515 | if (zmin > d->GetCoord1()) zmin=d->GetCoord1(); | |
516 | if (zmax < d->GetCoord1()) zmax=d->GetCoord1(); | |
e43c066c | 517 | // MI addition - find all labels in cluster |
518 | for (Int_t dlab=0;dlab<10;dlab++){ | |
babd135a | 519 | Int_t digitlab = (d->GetTracks())[dlab]; |
520 | if (digitlab<0) continue; | |
e43c066c | 521 | AddLabel(milab,digitlab); |
babd135a | 522 | } |
e43c066c | 523 | if (milab[9]>0) CheckLabels2(milab); |
524 | } | |
babd135a | 525 | CheckLabels2(milab); |
e43c066c | 526 | // |
527 | //Int_t idy = (fNlayer[fI]==0)? 2:3; | |
528 | //for (Int_t iz=zmin; iz<=zmax;iz+=2) | |
529 | //Int_t idy = (ymax-ymin)/4.; // max 2 clusters | |
530 | Int_t idy = 0; // max 2 clusters | |
531 | if (fNlayer[fI]==0 &&idy<3) idy=3; | |
532 | if (fNlayer[fI]==1 &&idy<4) idy=4; | |
533 | Int_t idz =3; | |
534 | for (Int_t iz=zmin; iz<=zmax;iz+=idz) | |
535 | for (Int_t iy=ymin; iy<=ymax;iy+=idy){ | |
536 | // | |
537 | Int_t ndigits =0; | |
0e986ef9 | 538 | Float_t y=0.,z=0.,q=0.; |
e43c066c | 539 | for (Int_t l=0; l<ni; l++) { |
540 | d=(AliITSdigitSPD*)digits->UncheckedAt(idx[l]); | |
541 | if (zmax-zmin>=idz || ymax-ymin>=idy){ | |
542 | if (TMath::Abs( d->GetCoord2()-iy)>0.75*idy) continue; | |
543 | if (TMath::Abs( d->GetCoord1()-iz)>0.75*idz) continue; | |
544 | } | |
545 | ndigits++; | |
0e986ef9 | 546 | Float_t qq=d->GetSignal(); |
e43c066c | 547 | y+=qq*fYSPD[d->GetCoord2()]; z+=qq*fZSPD[d->GetCoord1()]; q+=qq; |
548 | ||
549 | } | |
550 | if (ndigits==0) continue; | |
551 | y/=q; z/=q; | |
552 | y-=fHwSPD; z-=fHlSPD; | |
553 | ||
554 | Float_t lp[5]; | |
555 | lp[0]=-(-y+fYshift[fI]); if (fI<=fLastSPD1) lp[0]=-lp[0]; | |
556 | lp[1]= -z+fZshift[fI]; | |
557 | // Float_t factor=TMath::Max(double(ni-3.),1.5); | |
558 | Float_t factory=TMath::Max(ymax-ymin,1); | |
559 | Float_t factorz=TMath::Max(zmax-zmin,1); | |
560 | factory*= factory; | |
561 | factorz*= factorz; | |
562 | //lp[2]= (fYpitchSPD*fYpitchSPD/12.)*factory; | |
563 | //lp[3]= (fZ1pitchSPD*fZ1pitchSPD/12.)*factorz; | |
564 | lp[2]= (fYpitchSPD*fYpitchSPD/12.); | |
565 | lp[3]= (fZ1pitchSPD*fZ1pitchSPD/12.); | |
566 | //lp[4]= q; | |
567 | lp[4]= (zmax-zmin+1)*100 + (ymax-ymin+1); | |
568 | ||
569 | milab[3]=fNdet[fI]; | |
570 | d=(AliITSdigitSPD*)digits->UncheckedAt(idx[0]); | |
571 | Int_t info[3] = {ymax-ymin+1,zmax-zmin+1,fNlayer[fI]}; | |
572 | new (cl[n]) AliITSclusterV2(milab,lp,info); n++; | |
573 | } | |
7f4044f0 | 574 | } |
e43c066c | 575 | |
38bcdcc1 | 576 | delete [] bins; |
7f4044f0 | 577 | } |
578 | ||
c391f9d9 | 579 | void AliITSclustererV2::FindClustersSPD(AliITSRawStream* input, |
580 | TClonesArray** clusters) | |
581 | { | |
582 | //------------------------------------------------------------ | |
583 | // Actual SPD cluster finder for raw data | |
584 | //------------------------------------------------------------ | |
585 | ||
586 | Int_t nClustersSPD = 0; | |
587 | Int_t kNzBins = fNzSPD + 2; | |
588 | Int_t kNyBins = fNySPD + 2; | |
589 | Int_t kMaxBin = kNzBins * kNyBins; | |
5102d526 | 590 | AliBin *binsSPD = new AliBin[kMaxBin]; |
591 | AliBin *binsSPDInit = new AliBin[kMaxBin]; | |
592 | AliBin *bins = NULL; | |
c391f9d9 | 593 | |
594 | // read raw data input stream | |
595 | while (kTRUE) { | |
596 | Bool_t next = input->Next(); | |
597 | if (!next || input->IsNewModule()) { | |
598 | Int_t iModule = input->GetPrevModuleID(); | |
599 | ||
600 | // when all data from a module was read, search for clusters | |
601 | if (bins) { | |
602 | clusters[iModule] = new TClonesArray("AliITSclusterV2"); | |
603 | Int_t nClusters = 0; | |
604 | ||
605 | for (Int_t iBin = 0; iBin < kMaxBin; iBin++) { | |
606 | if (bins[iBin].IsUsed()) continue; | |
607 | Int_t nBins = 0; | |
608 | Int_t idxBins[200]; | |
609 | FindCluster(iBin, kNzBins, bins, nBins, idxBins); | |
610 | if (nBins == 200) { | |
611 | Error("FindClustersSPD", "SPD: Too big cluster !\n"); | |
612 | continue; | |
613 | } | |
614 | ||
615 | Int_t label[4]; | |
616 | label[0] = -2; | |
617 | label[1] = -2; | |
618 | label[2] = -2; | |
619 | // label[3] = iModule; | |
620 | label[3] = fNdet[iModule]; | |
621 | ||
622 | Int_t ymin = (idxBins[0] / kNzBins) - 1; | |
623 | Int_t ymax = ymin; | |
624 | Int_t zmin = (idxBins[0] % kNzBins) - 1; | |
625 | Int_t zmax = zmin; | |
c391f9d9 | 626 | for (Int_t idx = 0; idx < nBins; idx++) { |
627 | Int_t iy = (idxBins[idx] / kNzBins) - 1; | |
628 | Int_t iz = (idxBins[idx] % kNzBins) - 1; | |
629 | if (ymin > iy) ymin = iy; | |
630 | if (ymax < iy) ymax = iy; | |
631 | if (zmin > iz) zmin = iz; | |
632 | if (zmax < iz) zmax = iz; | |
c391f9d9 | 633 | } |
0e986ef9 | 634 | |
635 | Int_t idy = 0; // max 2 clusters | |
636 | if ((iModule <= fLastSPD1) &&idy<3) idy=3; | |
637 | if ((iModule > fLastSPD1) &&idy<4) idy=4; | |
638 | Int_t idz =3; | |
639 | for (Int_t iiz=zmin; iiz<=zmax;iiz+=idz) | |
640 | for (Int_t iiy=ymin; iiy<=ymax;iiy+=idy){ | |
641 | // | |
642 | Int_t ndigits =0; | |
643 | Float_t y=0.,z=0.,q=0.; | |
644 | for (Int_t idx = 0; idx < nBins; idx++) { | |
645 | Int_t iy = (idxBins[idx] / kNzBins) - 1; | |
646 | Int_t iz = (idxBins[idx] % kNzBins) - 1; | |
647 | if (zmax-zmin>=idz || ymax-ymin>=idy){ | |
648 | if (TMath::Abs(iy-iiy)>0.75*idy) continue; | |
649 | if (TMath::Abs(iz-iiz)>0.75*idz) continue; | |
650 | } | |
651 | ndigits++; | |
652 | Float_t qBin = bins[idxBins[idx]].GetQ(); | |
653 | y += qBin * fYSPD[iy]; | |
654 | z += qBin * fZSPD[iz]; | |
655 | q += qBin; | |
656 | } | |
657 | if (ndigits==0) continue; | |
658 | y /= q; | |
659 | z /= q; | |
660 | y -= fHwSPD; | |
661 | z -= fHlSPD; | |
662 | ||
663 | Float_t hit[5]; // y, z, sigma(y)^2, sigma(z)^2, charge | |
664 | hit[0] = -(-y+fYshift[iModule]); | |
665 | if (iModule <= fLastSPD1) hit[0] = -hit[0]; | |
666 | hit[1] = -z+fZshift[iModule]; | |
667 | hit[2] = fYpitchSPD*fYpitchSPD/12.; | |
668 | hit[3] = fZ1pitchSPD*fZ1pitchSPD/12.; | |
669 | // hit[4] = q; | |
670 | hit[4] = (zmax-zmin+1)*100 + (ymax-ymin+1); | |
671 | // CheckLabels(label); | |
672 | Int_t info[3]={ymax-ymin+1,zmax-zmin+1,fNlayer[iModule]}; | |
673 | new (clusters[iModule]->AddrAt(nClusters)) | |
674 | AliITSclusterV2(label, hit,info); | |
675 | nClusters++; | |
676 | } | |
c391f9d9 | 677 | } |
678 | ||
679 | nClustersSPD += nClusters; | |
5102d526 | 680 | bins = NULL; |
c391f9d9 | 681 | } |
682 | ||
683 | if (!next) break; | |
5102d526 | 684 | bins = binsSPD; |
685 | memcpy(binsSPD,binsSPDInit,sizeof(AliBin)*kMaxBin); | |
c391f9d9 | 686 | } |
687 | ||
688 | // fill the current digit into the bins array | |
689 | Int_t index = (input->GetCoord2()+1) * kNzBins + (input->GetCoord1()+1); | |
690 | bins[index].SetIndex(index); | |
691 | bins[index].SetMask(1); | |
692 | bins[index].SetQ(1); | |
693 | } | |
694 | ||
5102d526 | 695 | delete [] binsSPDInit; |
696 | delete [] binsSPD; | |
697 | ||
c391f9d9 | 698 | Info("FindClustersSPD", "found clusters in ITS SPD: %d", nClustersSPD); |
699 | } | |
700 | ||
701 | ||
7f4044f0 | 702 | Bool_t AliITSclustererV2::IsMaximum(Int_t k,Int_t max,const AliBin *bins) { |
703 | //------------------------------------------------------------ | |
704 | //is this a local maximum ? | |
705 | //------------------------------------------------------------ | |
706 | UShort_t q=bins[k].GetQ(); | |
707 | if (q==1023) return kFALSE; | |
708 | if (bins[k-max].GetQ() > q) return kFALSE; | |
709 | if (bins[k-1 ].GetQ() > q) return kFALSE; | |
710 | if (bins[k+max].GetQ() > q) return kFALSE; | |
711 | if (bins[k+1 ].GetQ() > q) return kFALSE; | |
712 | if (bins[k-max-1].GetQ() > q) return kFALSE; | |
713 | if (bins[k+max-1].GetQ() > q) return kFALSE; | |
714 | if (bins[k+max+1].GetQ() > q) return kFALSE; | |
715 | if (bins[k-max+1].GetQ() > q) return kFALSE; | |
716 | return kTRUE; | |
717 | } | |
718 | ||
719 | void AliITSclustererV2:: | |
720 | FindPeaks(Int_t k,Int_t max,AliBin *b,Int_t *idx,UInt_t *msk,Int_t& n) { | |
721 | //------------------------------------------------------------ | |
722 | //find local maxima | |
723 | //------------------------------------------------------------ | |
724 | if (n<31) | |
725 | if (IsMaximum(k,max,b)) { | |
726 | idx[n]=k; msk[n]=(2<<n); | |
727 | n++; | |
728 | } | |
729 | b[k].SetMask(0); | |
730 | if (b[k-max].GetMask()&1) FindPeaks(k-max,max,b,idx,msk,n); | |
731 | if (b[k-1 ].GetMask()&1) FindPeaks(k-1 ,max,b,idx,msk,n); | |
732 | if (b[k+max].GetMask()&1) FindPeaks(k+max,max,b,idx,msk,n); | |
733 | if (b[k+1 ].GetMask()&1) FindPeaks(k+1 ,max,b,idx,msk,n); | |
734 | } | |
735 | ||
736 | void AliITSclustererV2:: | |
737 | MarkPeak(Int_t k, Int_t max, AliBin *bins, UInt_t m) { | |
738 | //------------------------------------------------------------ | |
739 | //mark this peak | |
740 | //------------------------------------------------------------ | |
741 | UShort_t q=bins[k].GetQ(); | |
742 | ||
743 | bins[k].SetMask(bins[k].GetMask()|m); | |
744 | ||
745 | if (bins[k-max].GetQ() <= q) | |
746 | if ((bins[k-max].GetMask()&m) == 0) MarkPeak(k-max,max,bins,m); | |
747 | if (bins[k-1 ].GetQ() <= q) | |
748 | if ((bins[k-1 ].GetMask()&m) == 0) MarkPeak(k-1 ,max,bins,m); | |
749 | if (bins[k+max].GetQ() <= q) | |
750 | if ((bins[k+max].GetMask()&m) == 0) MarkPeak(k+max,max,bins,m); | |
751 | if (bins[k+1 ].GetQ() <= q) | |
752 | if ((bins[k+1 ].GetMask()&m) == 0) MarkPeak(k+1 ,max,bins,m); | |
753 | } | |
754 | ||
755 | void AliITSclustererV2:: | |
756 | MakeCluster(Int_t k,Int_t max,AliBin *bins,UInt_t m,AliITSclusterV2 &c) { | |
757 | //------------------------------------------------------------ | |
758 | //make cluster using digits of this peak | |
759 | //------------------------------------------------------------ | |
760 | Float_t q=(Float_t)bins[k].GetQ(); | |
761 | Int_t i=k/max, j=k-i*max; | |
762 | ||
763 | c.SetQ(c.GetQ()+q); | |
764 | c.SetY(c.GetY()+i*q); | |
765 | c.SetZ(c.GetZ()+j*q); | |
766 | c.SetSigmaY2(c.GetSigmaY2()+i*i*q); | |
767 | c.SetSigmaZ2(c.GetSigmaZ2()+j*j*q); | |
768 | ||
769 | bins[k].SetMask(0xFFFFFFFE); | |
770 | ||
771 | if (bins[k-max].GetMask() == m) MakeCluster(k-max,max,bins,m,c); | |
772 | if (bins[k-1 ].GetMask() == m) MakeCluster(k-1 ,max,bins,m,c); | |
773 | if (bins[k+max].GetMask() == m) MakeCluster(k+max,max,bins,m,c); | |
774 | if (bins[k+1 ].GetMask() == m) MakeCluster(k+1 ,max,bins,m,c); | |
775 | } | |
776 | ||
777 | void AliITSclustererV2:: | |
c391f9d9 | 778 | FindClustersSDD(AliBin* bins[2], Int_t nMaxBin, Int_t nzBins, |
779 | const TClonesArray *digits, TClonesArray *clusters) { | |
7f4044f0 | 780 | //------------------------------------------------------------ |
781 | // Actual SDD cluster finder | |
782 | //------------------------------------------------------------ | |
7f4044f0 | 783 | Int_t ncl=0; TClonesArray &cl=*clusters; |
784 | for (Int_t s=0; s<2; s++) | |
c391f9d9 | 785 | for (Int_t i=0; i<nMaxBin; i++) { |
7f4044f0 | 786 | if (bins[s][i].IsUsed()) continue; |
787 | Int_t idx[32]; UInt_t msk[32]; Int_t npeaks=0; | |
c391f9d9 | 788 | FindPeaks(i, nzBins, bins[s], idx, msk, npeaks); |
7f4044f0 | 789 | |
790 | if (npeaks>30) continue; | |
e43c066c | 791 | if (npeaks==0) continue; |
7f4044f0 | 792 | |
793 | Int_t k,l; | |
794 | for (k=0; k<npeaks-1; k++){//mark adjacent peaks | |
795 | if (idx[k] < 0) continue; //this peak is already removed | |
796 | for (l=k+1; l<npeaks; l++) { | |
797 | if (idx[l] < 0) continue; //this peak is already removed | |
c391f9d9 | 798 | Int_t ki=idx[k]/nzBins, kj=idx[k] - ki*nzBins; |
799 | Int_t li=idx[l]/nzBins, lj=idx[l] - li*nzBins; | |
7f4044f0 | 800 | Int_t di=TMath::Abs(ki - li); |
801 | Int_t dj=TMath::Abs(kj - lj); | |
802 | if (di>1 || dj>1) continue; | |
803 | if (bins[s][idx[k]].GetQ() > bins[s][idx[l]].GetQ()) { | |
804 | msk[l]=msk[k]; | |
805 | idx[l]*=-1; | |
806 | } else { | |
807 | msk[k]=msk[l]; | |
808 | idx[k]*=-1; | |
809 | break; | |
810 | } | |
811 | } | |
812 | } | |
813 | ||
814 | for (k=0; k<npeaks; k++) { | |
c391f9d9 | 815 | MarkPeak(TMath::Abs(idx[k]), nzBins, bins[s], msk[k]); |
7f4044f0 | 816 | } |
817 | ||
818 | for (k=0; k<npeaks; k++) { | |
819 | if (idx[k] < 0) continue; //removed peak | |
820 | AliITSclusterV2 c; | |
c391f9d9 | 821 | MakeCluster(idx[k], nzBins, bins[s], msk[k], c); |
babd135a | 822 | //mi change |
823 | Int_t milab[10]; | |
824 | for (Int_t ilab=0;ilab<10;ilab++){ | |
825 | milab[ilab]=-2; | |
826 | } | |
e43c066c | 827 | Int_t maxi=0,mini=0,maxj=0,minj=0; |
828 | //AliBin *bmax=&bins[s][idx[k]]; | |
829 | //Float_t max = TMath::Max(TMath::Abs(bmax->GetQ())/5.,3.); | |
830 | Float_t max=3; | |
831 | for (Int_t di=-2; di<=2;di++) | |
832 | for (Int_t dj=-3;dj<=3;dj++){ | |
833 | Int_t index = idx[k]+di+dj*nzBins; | |
834 | if (index<0) continue; | |
835 | if (index>=nMaxBin) continue; | |
836 | AliBin *b=&bins[s][index]; | |
837 | if (TMath::Abs(b->GetQ())>max){ | |
838 | if (di>maxi) maxi=di; | |
839 | if (di<mini) mini=di; | |
840 | if (dj>maxj) maxj=dj; | |
841 | if (dj<minj) minj=dj; | |
842 | // | |
0e986ef9 | 843 | if(digits) { |
844 | if (TMath::Abs(di)<2&&TMath::Abs(dj)<2){ | |
845 | AliITSdigitSDD* d=(AliITSdigitSDD*)digits->UncheckedAt(b->GetIndex()); | |
846 | for (Int_t itrack=0;itrack<10;itrack++){ | |
847 | Int_t track = (d->GetTracks())[itrack]; | |
848 | if (track>=0) { | |
849 | AddLabel(milab, track); | |
850 | } | |
e43c066c | 851 | } |
852 | } | |
853 | } | |
854 | } | |
855 | } | |
856 | ||
857 | /* | |
858 | Float_t s2 = c.GetSigmaY2()/c.GetQ() - c.GetY()*c.GetY(); | |
859 | Float_t w=par->GetPadPitchWidth(sec); | |
860 | c.SetSigmaY2(s2); | |
861 | if (s2 != 0.) { | |
862 | c.SetSigmaY2(c.GetSigmaY2()*0.108); | |
863 | if (sec<par->GetNInnerSector()) c.SetSigmaY2(c.GetSigmaY2()*2.07); | |
864 | } | |
865 | s2 = c.GetSigmaZ2()/c.GetQ() - c.GetZ()*c.GetZ(); | |
866 | w=par->GetZWidth(); | |
867 | c.SetSigmaZ2(s2); | |
868 | ||
869 | if (s2 != 0.) { | |
870 | c.SetSigmaZ2(c.GetSigmaZ2()*0.169); | |
871 | if (sec<par->GetNInnerSector()) c.SetSigmaZ2(c.GetSigmaZ2()*1.77); | |
872 | } | |
7f4044f0 | 873 | */ |
874 | ||
875 | c.SetSigmaY2(0.0030*0.0030); | |
876 | c.SetSigmaZ2(0.0020*0.0020); | |
877 | c.SetDetectorIndex(fNdet[fI]); | |
878 | ||
879 | Float_t y=c.GetY(),z=c.GetZ(), q=c.GetQ(); | |
880 | y/=q; z/=q; | |
e43c066c | 881 | // |
882 | //Float_t s2 = c.GetSigmaY2()/c.GetQ() - y*y; | |
883 | // c.SetSigmaY2(s2); | |
884 | //s2 = c.GetSigmaZ2()/c.GetQ() - z*z; | |
885 | //c.SetSigmaZ2(s2); | |
886 | // | |
7f4044f0 | 887 | y=(y-0.5)*fYpitchSDD; |
888 | y-=fHwSDD; | |
889 | y-=fYoffSDD; //delay ? | |
890 | if (s) y=-y; | |
891 | ||
892 | z=(z-0.5)*fZpitchSDD; | |
893 | z-=fHlSDD; | |
894 | ||
f363d377 | 895 | y=-(-y+fYshift[fI]); |
896 | z= -z+fZshift[fI]; | |
7f4044f0 | 897 | c.SetY(y); |
898 | c.SetZ(z); | |
e43c066c | 899 | c.SetNy(maxj-minj+1); |
900 | c.SetNz(maxi-mini+1); | |
901 | c.SetType(npeaks); | |
e3d91d99 | 902 | c.SetQ(q/12.7); //to be consistent with the SSD charges |
903 | ||
904 | if (c.GetQ() < 20.) continue; //noise cluster | |
e43c066c | 905 | |
906 | if (digits) { | |
907 | // AliBin *b=&bins[s][idx[k]]; | |
908 | // AliITSdigitSDD* d=(AliITSdigitSDD*)digits->UncheckedAt(b->GetIndex()); | |
c391f9d9 | 909 | { |
e43c066c | 910 | //Int_t lab[3]; |
911 | //lab[0]=(d->GetTracks())[0]; | |
912 | //lab[1]=(d->GetTracks())[1]; | |
913 | //lab[2]=(d->GetTracks())[2]; | |
914 | //CheckLabels(lab); | |
babd135a | 915 | CheckLabels2(milab); |
916 | c.SetLabel(milab[0],0); | |
917 | c.SetLabel(milab[1],1); | |
918 | c.SetLabel(milab[2],2); | |
e43c066c | 919 | c.SetLayer(fNlayer[fI]); |
c391f9d9 | 920 | } |
921 | } | |
7f4044f0 | 922 | new (cl[ncl]) AliITSclusterV2(c); ncl++; |
923 | } | |
924 | } | |
c391f9d9 | 925 | } |
926 | ||
927 | void AliITSclustererV2:: | |
928 | FindClustersSDD(const TClonesArray *digits, TClonesArray *clusters) { | |
929 | //------------------------------------------------------------ | |
930 | // Actual SDD cluster finder | |
931 | //------------------------------------------------------------ | |
932 | Int_t kNzBins = fNzSDD + 2; | |
933 | const Int_t kMAXBIN=kNzBins*(fNySDD+2); | |
934 | ||
935 | AliBin *bins[2]; | |
936 | bins[0]=new AliBin[kMAXBIN]; | |
937 | bins[1]=new AliBin[kMAXBIN]; | |
938 | ||
939 | AliITSdigitSDD *d=0; | |
940 | Int_t i, ndigits=digits->GetEntriesFast(); | |
941 | for (i=0; i<ndigits; i++) { | |
942 | d=(AliITSdigitSDD*)digits->UncheckedAt(i); | |
943 | Int_t y=d->GetCoord2()+1; //y | |
944 | Int_t z=d->GetCoord1()+1; //z | |
945 | Int_t q=d->GetSignal(); | |
e43c066c | 946 | if (q<3) continue; |
e3d91d99 | 947 | |
c391f9d9 | 948 | if (z <= fNzSDD) { |
949 | bins[0][y*kNzBins+z].SetQ(q); | |
950 | bins[0][y*kNzBins+z].SetMask(1); | |
951 | bins[0][y*kNzBins+z].SetIndex(i); | |
952 | } else { | |
953 | z-=fNzSDD; | |
954 | bins[1][y*kNzBins+z].SetQ(q); | |
955 | bins[1][y*kNzBins+z].SetMask(1); | |
956 | bins[1][y*kNzBins+z].SetIndex(i); | |
957 | } | |
958 | } | |
959 | ||
960 | FindClustersSDD(bins, kMAXBIN, kNzBins, digits, clusters); | |
7f4044f0 | 961 | |
962 | delete[] bins[0]; | |
963 | delete[] bins[1]; | |
0e986ef9 | 964 | |
7f4044f0 | 965 | } |
966 | ||
c391f9d9 | 967 | void AliITSclustererV2::FindClustersSDD(AliITSRawStream* input, |
968 | TClonesArray** clusters) | |
969 | { | |
970 | //------------------------------------------------------------ | |
971 | // Actual SDD cluster finder for raw data | |
972 | //------------------------------------------------------------ | |
973 | Int_t nClustersSDD = 0; | |
974 | Int_t kNzBins = fNzSDD + 2; | |
975 | Int_t kMaxBin = kNzBins * (fNySDD+2); | |
5102d526 | 976 | AliBin *binsSDDInit = new AliBin[kMaxBin]; |
977 | AliBin *binsSDD1 = new AliBin[kMaxBin]; | |
978 | AliBin *binsSDD2 = new AliBin[kMaxBin]; | |
979 | AliBin *bins[2] = {NULL, NULL}; | |
c391f9d9 | 980 | |
981 | // read raw data input stream | |
982 | while (kTRUE) { | |
983 | Bool_t next = input->Next(); | |
984 | if (!next || input->IsNewModule()) { | |
985 | Int_t iModule = input->GetPrevModuleID(); | |
986 | ||
987 | // when all data from a module was read, search for clusters | |
988 | if (bins[0]) { | |
989 | clusters[iModule] = new TClonesArray("AliITSclusterV2"); | |
990 | fI = iModule; | |
991 | FindClustersSDD(bins, kMaxBin, kNzBins, NULL, clusters[iModule]); | |
992 | Int_t nClusters = clusters[iModule]->GetEntriesFast(); | |
993 | nClustersSDD += nClusters; | |
5102d526 | 994 | bins[0] = bins[1] = NULL; |
c391f9d9 | 995 | } |
996 | ||
997 | if (!next) break; | |
5102d526 | 998 | bins[0]=binsSDD1; |
999 | bins[1]=binsSDD2; | |
1000 | memcpy(binsSDD1,binsSDDInit,sizeof(AliBin)*kMaxBin); | |
1001 | memcpy(binsSDD2,binsSDDInit,sizeof(AliBin)*kMaxBin); | |
c391f9d9 | 1002 | } |
1003 | ||
1004 | // fill the current digit into the bins array | |
0e986ef9 | 1005 | if(input->GetSignal()>=3) { |
1006 | Int_t iz = input->GetCoord1()+1; | |
1007 | Int_t side = ((iz <= fNzSDD) ? 0 : 1); | |
1008 | iz -= side*fNzSDD; | |
1009 | Int_t index = (input->GetCoord2()+1) * kNzBins + iz; | |
1010 | bins[side][index].SetQ(input->GetSignal()); | |
1011 | bins[side][index].SetMask(1); | |
1012 | bins[side][index].SetIndex(index); | |
1013 | } | |
c391f9d9 | 1014 | } |
1015 | ||
5102d526 | 1016 | delete[] binsSDD1; |
1017 | delete[] binsSDD2; | |
1018 | delete[] binsSDDInit; | |
1019 | ||
c391f9d9 | 1020 | Info("FindClustersSDD", "found clusters in ITS SDD: %d", nClustersSDD); |
1021 | } | |
1022 | ||
e43c066c | 1023 | |
1024 | ||
c391f9d9 | 1025 | void AliITSclustererV2:: |
1026 | FindClustersSSD(Ali1Dcluster* neg, Int_t nn, | |
1027 | Ali1Dcluster* pos, Int_t np, | |
1028 | TClonesArray *clusters) { | |
1029 | //------------------------------------------------------------ | |
1030 | // Actual SSD cluster finder | |
1031 | //------------------------------------------------------------ | |
1032 | TClonesArray &cl=*clusters; | |
e43c066c | 1033 | // |
c391f9d9 | 1034 | Float_t tanp=fTanP, tann=fTanN; |
1035 | if (fI>fLastSSD1) {tann=fTanP; tanp=fTanN;} | |
c391f9d9 | 1036 | Int_t idet=fNdet[fI]; |
1037 | Int_t ncl=0; | |
e43c066c | 1038 | // |
1039 | Int_t negativepair[30000]; | |
1040 | Int_t cnegative[3000]; | |
1041 | Int_t cused1[3000]; | |
1042 | Int_t positivepair[30000]; | |
1043 | Int_t cpositive[3000]; | |
1044 | Int_t cused2[3000]; | |
1045 | for (Int_t i=0;i<3000;i++) {cnegative[i]=0; cused1[i]=0;} | |
1046 | for (Int_t i=0;i<3000;i++) {cpositive[i]=0; cused2[i]=0;} | |
0774985e | 1047 | static Short_t pairs[1000][1000]; |
1048 | memset(pairs,0,sizeof(Short_t)*1000000); | |
1049 | // Short_t ** pairs = new Short_t*[1000]; | |
1050 | // for (Int_t i=0; i<1000; i++) { | |
1051 | // pairs[i] = new Short_t[1000]; | |
1052 | // memset(pairs[i],0,sizeof(Short_t)*1000); | |
1053 | // } | |
e43c066c | 1054 | // |
1055 | // find available pairs | |
1056 | // | |
1057 | for (Int_t i=0; i<np; i++) { | |
1058 | Float_t yp=pos[i].GetY()*fYpitchSSD; | |
1059 | if (pos[i].GetQ()<3) continue; | |
1060 | for (Int_t j=0; j<nn; j++) { | |
1061 | if (neg[j].GetQ()<3) continue; | |
1062 | Float_t yn=neg[j].GetY()*fYpitchSSD; | |
1063 | Float_t zt=(2*fHlSSD*tanp + yp - yn)/(tann+tanp); | |
1064 | Float_t yt=yn + tann*zt; | |
1065 | zt-=fHlSSD; yt-=fHwSSD; | |
1066 | if (TMath::Abs(yt)<fHwSSD+0.01) | |
1067 | if (TMath::Abs(zt)<fHlSSD+0.01*(neg[j].GetNd()+pos[i].GetNd())) { | |
1068 | negativepair[i*10+cnegative[i]] =j; //index | |
1069 | positivepair[j*10+cpositive[j]] =i; | |
1070 | cnegative[i]++; //counters | |
1071 | cpositive[j]++; | |
1072 | pairs[i][j]=100; | |
1073 | } | |
1074 | } | |
1075 | } | |
1076 | // | |
1077 | for (Int_t i=0; i<np; i++) { | |
1078 | Float_t yp=pos[i].GetY()*fYpitchSSD; | |
1079 | if (pos[i].GetQ()<3) continue; | |
1080 | for (Int_t j=0; j<nn; j++) { | |
1081 | if (neg[j].GetQ()<3) continue; | |
1082 | if (cpositive[j]&&cnegative[i]) continue; | |
1083 | Float_t yn=neg[j].GetY()*fYpitchSSD; | |
1084 | Float_t zt=(2*fHlSSD*tanp + yp - yn)/(tann+tanp); | |
1085 | Float_t yt=yn + tann*zt; | |
1086 | zt-=fHlSSD; yt-=fHwSSD; | |
1087 | if (TMath::Abs(yt)<fHwSSD+0.1) | |
1088 | if (TMath::Abs(zt)<fHlSSD+0.15) { | |
1089 | if (cnegative[i]==0) pos[i].SetNd(100); // not available pair | |
1090 | if (cpositive[j]==0) neg[j].SetNd(100); // not available pair | |
1091 | negativepair[i*10+cnegative[i]] =j; //index | |
1092 | positivepair[j*10+cpositive[j]] =i; | |
1093 | cnegative[i]++; //counters | |
1094 | cpositive[j]++; | |
1095 | pairs[i][j]=100; | |
1096 | } | |
1097 | } | |
1098 | } | |
1099 | // | |
1100 | Float_t lp[5]; | |
1101 | Int_t milab[10]; | |
1102 | Double_t ratio; | |
1103 | ||
1104 | // | |
1105 | // sign gold tracks | |
1106 | // | |
1107 | for (Int_t ip=0;ip<np;ip++){ | |
1108 | Float_t ybest=1000,zbest=1000,qbest=0; | |
1109 | // | |
1110 | // select gold clusters | |
1111 | if ( (cnegative[ip]==1) && cpositive[negativepair[10*ip]]==1){ | |
1112 | Float_t yp=pos[ip].GetY()*fYpitchSSD; | |
1113 | Int_t j = negativepair[10*ip]; | |
1114 | ratio = (pos[ip].GetQ()-neg[j].GetQ())/(pos[ip].GetQ()+neg[j].GetQ()); | |
1115 | // | |
1116 | Float_t yn=neg[j].GetY()*fYpitchSSD; | |
1117 | Float_t zt=(2*fHlSSD*tanp + yp - yn)/(tann+tanp); | |
1118 | Float_t yt=yn + tann*zt; | |
1119 | zt-=fHlSSD; yt-=fHwSSD; | |
1120 | ybest=yt; zbest=zt; | |
1121 | qbest=0.5*(pos[ip].GetQ()+neg[j].GetQ()); | |
1122 | lp[0]=-(-ybest+fYshift[fI]); | |
1123 | lp[1]= -zbest+fZshift[fI]; | |
1124 | lp[2]=0.0025*0.0025; //SigmaY2 | |
1125 | lp[3]=0.110*0.110; //SigmaZ2 | |
1126 | ||
1127 | lp[4]=qbest; //Q | |
1128 | for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2; | |
1129 | for (Int_t ilab=0;ilab<3;ilab++){ | |
1130 | milab[ilab] = pos[ip].GetLabel(ilab); | |
1131 | milab[ilab+3] = neg[j].GetLabel(ilab); | |
1132 | } | |
1133 | // | |
1134 | CheckLabels2(milab); | |
1135 | milab[3]=(((ip<<10) + j)<<10) + idet; // pos|neg|det | |
1136 | Int_t info[3] = {pos[ip].GetNd(),neg[j].GetNd(),fNlayer[fI]}; | |
1137 | AliITSclusterV2 * cl2 = new (cl[ncl]) AliITSclusterV2(milab,lp,info); | |
1138 | ncl++; | |
1139 | cl2->SetChargeRatio(ratio); | |
1140 | cl2->SetType(1); | |
1141 | pairs[ip][j]=1; | |
1142 | if ((pos[ip].GetNd()+neg[j].GetNd())>6){ //multi cluster | |
1143 | cl2->SetType(2); | |
1144 | pairs[ip][j]=2; | |
1145 | } | |
1146 | cused1[ip]++; | |
1147 | cused2[j]++; | |
1148 | } | |
1149 | } | |
1150 | ||
1151 | for (Int_t ip=0;ip<np;ip++){ | |
1152 | Float_t ybest=1000,zbest=1000,qbest=0; | |
1153 | // | |
1154 | // | |
1155 | // select "silber" cluster | |
1156 | if ( cnegative[ip]==1 && cpositive[negativepair[10*ip]]==2){ | |
1157 | Int_t in = negativepair[10*ip]; | |
1158 | Int_t ip2 = positivepair[10*in]; | |
1159 | if (ip2==ip) ip2 = positivepair[10*in+1]; | |
1160 | Float_t pcharge = pos[ip].GetQ()+pos[ip2].GetQ(); | |
1161 | if (TMath::Abs(pcharge-neg[in].GetQ())<10){ | |
1162 | // | |
1163 | // add first pair | |
1164 | if (pairs[ip][in]==100){ // | |
1165 | Float_t yp=pos[ip].GetY()*fYpitchSSD; | |
1166 | Float_t yn=neg[in].GetY()*fYpitchSSD; | |
1167 | Float_t zt=(2*fHlSSD*tanp + yp - yn)/(tann+tanp); | |
1168 | Float_t yt=yn + tann*zt; | |
1169 | zt-=fHlSSD; yt-=fHwSSD; | |
1170 | ybest =yt; zbest=zt; | |
1171 | qbest =pos[ip].GetQ(); | |
1172 | lp[0]=-(-ybest+fYshift[fI]); | |
1173 | lp[1]= -zbest+fZshift[fI]; | |
1174 | lp[2]=0.0025*0.0025; //SigmaY2 | |
1175 | lp[3]=0.110*0.110; //SigmaZ2 | |
1176 | ||
1177 | lp[4]=qbest; //Q | |
1178 | for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2; | |
1179 | for (Int_t ilab=0;ilab<3;ilab++){ | |
1180 | milab[ilab] = pos[ip].GetLabel(ilab); | |
1181 | milab[ilab+3] = neg[in].GetLabel(ilab); | |
1182 | } | |
1183 | // | |
1184 | CheckLabels2(milab); | |
1185 | ratio = (pos[ip].GetQ()-neg[in].GetQ())/(pos[ip].GetQ()+neg[in].GetQ()); | |
1186 | milab[3]=(((ip<<10) + in)<<10) + idet; // pos|neg|det | |
1187 | Int_t info[3] = {pos[ip].GetNd(),neg[in].GetNd(),fNlayer[fI]}; | |
1188 | AliITSclusterV2 * cl2 = new (cl[ncl]) AliITSclusterV2(milab,lp,info); | |
1189 | ncl++; | |
1190 | cl2->SetChargeRatio(ratio); | |
1191 | cl2->SetType(5); | |
1192 | pairs[ip][in] = 5; | |
1193 | if ((pos[ip].GetNd()+neg[in].GetNd())>6){ //multi cluster | |
1194 | cl2->SetType(6); | |
1195 | pairs[ip][in] = 6; | |
1196 | } | |
1197 | } | |
1198 | // | |
1199 | // add second pair | |
1200 | ||
1201 | // if (!(cused1[ip2] || cused2[in])){ // | |
1202 | if (pairs[ip2][in]==100){ | |
1203 | Float_t yp=pos[ip2].GetY()*fYpitchSSD; | |
1204 | Float_t yn=neg[in].GetY()*fYpitchSSD; | |
1205 | Float_t zt=(2*fHlSSD*tanp + yp - yn)/(tann+tanp); | |
1206 | Float_t yt=yn + tann*zt; | |
1207 | zt-=fHlSSD; yt-=fHwSSD; | |
1208 | ybest =yt; zbest=zt; | |
1209 | qbest =pos[ip2].GetQ(); | |
1210 | lp[0]=-(-ybest+fYshift[fI]); | |
1211 | lp[1]= -zbest+fZshift[fI]; | |
1212 | lp[2]=0.0025*0.0025; //SigmaY2 | |
1213 | lp[3]=0.110*0.110; //SigmaZ2 | |
1214 | ||
1215 | lp[4]=qbest; //Q | |
1216 | for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2; | |
1217 | for (Int_t ilab=0;ilab<3;ilab++){ | |
1218 | milab[ilab] = pos[ip2].GetLabel(ilab); | |
1219 | milab[ilab+3] = neg[in].GetLabel(ilab); | |
1220 | } | |
1221 | // | |
1222 | CheckLabels2(milab); | |
1223 | ratio = (pos[ip2].GetQ()-neg[in].GetQ())/(pos[ip2].GetQ()+neg[in].GetQ()); | |
1224 | milab[3]=(((ip2<<10) + in)<<10) + idet; // pos|neg|det | |
1225 | Int_t info[3] = {pos[ip2].GetNd(),neg[in].GetNd(),fNlayer[fI]}; | |
1226 | AliITSclusterV2 *cl2 = new (cl[ncl]) AliITSclusterV2(milab,lp,info); | |
1227 | ncl++; | |
1228 | cl2->SetChargeRatio(ratio); | |
1229 | cl2->SetType(5); | |
1230 | pairs[ip2][in] =5; | |
1231 | if ((pos[ip2].GetNd()+neg[in].GetNd())>6){ //multi cluster | |
1232 | cl2->SetType(6); | |
1233 | pairs[ip2][in] =6; | |
1234 | } | |
1235 | } | |
1236 | cused1[ip]++; | |
1237 | cused1[ip2]++; | |
1238 | cused2[in]++; | |
1239 | } | |
1240 | } | |
1241 | } | |
1242 | ||
1243 | // | |
1244 | for (Int_t jn=0;jn<nn;jn++){ | |
1245 | if (cused2[jn]) continue; | |
1246 | Float_t ybest=1000,zbest=1000,qbest=0; | |
1247 | // select "silber" cluster | |
1248 | if ( cpositive[jn]==1 && cnegative[positivepair[10*jn]]==2){ | |
1249 | Int_t ip = positivepair[10*jn]; | |
1250 | Int_t jn2 = negativepair[10*ip]; | |
1251 | if (jn2==jn) jn2 = negativepair[10*ip+1]; | |
1252 | Float_t pcharge = neg[jn].GetQ()+neg[jn2].GetQ(); | |
1253 | // | |
1254 | if (TMath::Abs(pcharge-pos[ip].GetQ())<10){ | |
1255 | // | |
1256 | // add first pair | |
1257 | // if (!(cused1[ip]||cused2[jn])){ | |
1258 | if (pairs[ip][jn]==100){ | |
1259 | Float_t yn=neg[jn].GetY()*fYpitchSSD; | |
1260 | Float_t yp=pos[ip].GetY()*fYpitchSSD; | |
1261 | Float_t zt=(2*fHlSSD*tanp + yp - yn)/(tann+tanp); | |
1262 | Float_t yt=yn + tann*zt; | |
1263 | zt-=fHlSSD; yt-=fHwSSD; | |
1264 | ybest =yt; zbest=zt; | |
1265 | qbest =neg[jn].GetQ(); | |
1266 | lp[0]=-(-ybest+fYshift[fI]); | |
1267 | lp[1]= -zbest+fZshift[fI]; | |
1268 | lp[2]=0.0025*0.0025; //SigmaY2 | |
1269 | lp[3]=0.110*0.110; //SigmaZ2 | |
1270 | ||
1271 | lp[4]=qbest; //Q | |
1272 | for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2; | |
1273 | for (Int_t ilab=0;ilab<3;ilab++){ | |
1274 | milab[ilab] = pos[ip].GetLabel(ilab); | |
1275 | milab[ilab+3] = neg[jn].GetLabel(ilab); | |
1276 | } | |
1277 | // | |
1278 | CheckLabels2(milab); | |
1279 | ratio = (pos[ip].GetQ()-neg[jn].GetQ())/(pos[ip].GetQ()+neg[jn].GetQ()); | |
1280 | milab[3]=(((ip<<10) + jn)<<10) + idet; // pos|neg|det | |
1281 | Int_t info[3] = {pos[ip].GetNd(),neg[jn].GetNd(),fNlayer[fI]}; | |
1282 | AliITSclusterV2 * cl2 = new (cl[ncl]) AliITSclusterV2(milab,lp,info); | |
1283 | ncl++; | |
1284 | cl2->SetChargeRatio(ratio); | |
1285 | cl2->SetType(7); | |
1286 | pairs[ip][jn] =7; | |
1287 | if ((pos[ip].GetNd()+neg[jn].GetNd())>6){ //multi cluster | |
1288 | cl2->SetType(8); | |
1289 | pairs[ip][jn]=8; | |
1290 | } | |
1291 | } | |
1292 | // | |
1293 | // add second pair | |
1294 | // if (!(cused1[ip]||cused2[jn2])){ | |
1295 | if (pairs[ip][jn2]==100){ | |
1296 | Float_t yn=neg[jn2].GetY()*fYpitchSSD; | |
1297 | Double_t yp=pos[ip].GetY()*fYpitchSSD; | |
1298 | Double_t zt=(2*fHlSSD*tanp + yp - yn)/(tann+tanp); | |
1299 | Double_t yt=yn + tann*zt; | |
1300 | zt-=fHlSSD; yt-=fHwSSD; | |
1301 | ybest =yt; zbest=zt; | |
1302 | qbest =neg[jn2].GetQ(); | |
1303 | lp[0]=-(-ybest+fYshift[fI]); | |
1304 | lp[1]= -zbest+fZshift[fI]; | |
1305 | lp[2]=0.0025*0.0025; //SigmaY2 | |
1306 | lp[3]=0.110*0.110; //SigmaZ2 | |
1307 | ||
1308 | lp[4]=qbest; //Q | |
1309 | for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2; | |
1310 | for (Int_t ilab=0;ilab<3;ilab++){ | |
1311 | milab[ilab] = pos[ip].GetLabel(ilab); | |
1312 | milab[ilab+3] = neg[jn2].GetLabel(ilab); | |
1313 | } | |
1314 | // | |
1315 | CheckLabels2(milab); | |
1316 | ratio = (pos[ip].GetQ()-neg[jn2].GetQ())/(pos[ip].GetQ()+neg[jn2].GetQ()); | |
1317 | milab[3]=(((ip<<10) + jn2)<<10) + idet; // pos|neg|det | |
1318 | Int_t info[3] = {pos[ip].GetNd(),neg[jn2].GetNd(),fNlayer[fI]}; | |
1319 | AliITSclusterV2* cl2 = new (cl[ncl]) AliITSclusterV2(milab,lp,info); | |
1320 | ncl++; | |
1321 | cl2->SetChargeRatio(ratio); | |
1322 | pairs[ip][jn2]=7; | |
1323 | cl2->SetType(7); | |
1324 | if ((pos[ip].GetNd()+neg[jn2].GetNd())>6){ //multi cluster | |
1325 | cl2->SetType(8); | |
1326 | pairs[ip][jn2]=8; | |
1327 | } | |
1328 | } | |
1329 | cused1[ip]++; | |
1330 | cused2[jn]++; | |
1331 | cused2[jn2]++; | |
1332 | } | |
1333 | } | |
1334 | } | |
1335 | ||
1336 | for (Int_t ip=0;ip<np;ip++){ | |
1337 | Float_t ybest=1000,zbest=1000,qbest=0; | |
1338 | // | |
1339 | // 2x2 clusters | |
1340 | // | |
1341 | if ( (cnegative[ip]<5) && cpositive[negativepair[10*ip]]<5){ | |
1342 | Float_t minchargediff =4.; | |
1343 | Int_t j=-1; | |
1344 | for (Int_t di=0;di<cnegative[ip];di++){ | |
1345 | Int_t jc = negativepair[ip*10+di]; | |
1346 | Float_t chargedif = pos[ip].GetQ()-neg[jc].GetQ(); | |
1347 | if (TMath::Abs(chargedif)<minchargediff){ | |
1348 | j =jc; | |
1349 | minchargediff = TMath::Abs(chargedif); | |
1350 | } | |
1351 | } | |
1352 | if (j<0) continue; // not proper cluster | |
1353 | Int_t count =0; | |
1354 | for (Int_t di=0;di<cnegative[ip];di++){ | |
1355 | Int_t jc = negativepair[ip*10+di]; | |
1356 | Float_t chargedif = pos[ip].GetQ()-neg[jc].GetQ(); | |
1357 | if (TMath::Abs(chargedif)<minchargediff+3.) count++; | |
1358 | } | |
1359 | if (count>1) continue; // more than one "proper" cluster for positive | |
1360 | // | |
1361 | count =0; | |
1362 | for (Int_t dj=0;dj<cpositive[j];dj++){ | |
1363 | Int_t ic = positivepair[j*10+dj]; | |
1364 | Float_t chargedif = pos[ic].GetQ()-neg[j].GetQ(); | |
1365 | if (TMath::Abs(chargedif)<minchargediff+3.) count++; | |
1366 | } | |
1367 | if (count>1) continue; // more than one "proper" cluster for negative | |
1368 | ||
1369 | Int_t jp = 0; | |
1370 | ||
1371 | count =0; | |
1372 | for (Int_t dj=0;dj<cnegative[jp];dj++){ | |
1373 | Int_t ic = positivepair[jp*10+dj]; | |
1374 | Float_t chargedif = pos[ic].GetQ()-neg[jp].GetQ(); | |
1375 | if (TMath::Abs(chargedif)<minchargediff+4.) count++; | |
1376 | } | |
1377 | if (count>1) continue; | |
1378 | if (pairs[ip][j]<100) continue; | |
1379 | // | |
1380 | //almost gold clusters | |
1381 | Float_t yp=pos[ip].GetY()*fYpitchSSD; | |
1382 | Float_t yn=neg[j].GetY()*fYpitchSSD; | |
1383 | Float_t zt=(2*fHlSSD*tanp + yp - yn)/(tann+tanp); | |
1384 | Float_t yt=yn + tann*zt; | |
1385 | zt-=fHlSSD; yt-=fHwSSD; | |
1386 | ybest=yt; zbest=zt; | |
1387 | qbest=0.5*(pos[ip].GetQ()+neg[j].GetQ()); | |
1388 | lp[0]=-(-ybest+fYshift[fI]); | |
1389 | lp[1]= -zbest+fZshift[fI]; | |
1390 | lp[2]=0.0025*0.0025; //SigmaY2 | |
1391 | lp[3]=0.110*0.110; //SigmaZ2 | |
1392 | lp[4]=qbest; //Q | |
1393 | for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2; | |
1394 | for (Int_t ilab=0;ilab<3;ilab++){ | |
1395 | milab[ilab] = pos[ip].GetLabel(ilab); | |
1396 | milab[ilab+3] = neg[j].GetLabel(ilab); | |
1397 | } | |
1398 | // | |
1399 | CheckLabels2(milab); | |
1400 | ratio = (pos[ip].GetQ()-neg[j].GetQ())/(pos[ip].GetQ()+neg[j].GetQ()); | |
1401 | milab[3]=(((ip<<10) + j)<<10) + idet; // pos|neg|det | |
1402 | Int_t info[3] = {pos[ip].GetNd(),neg[j].GetNd(),fNlayer[fI]}; | |
1403 | AliITSclusterV2 * cl2 = new (cl[ncl]) AliITSclusterV2(milab,lp,info); | |
1404 | ncl++; | |
1405 | cl2->SetChargeRatio(ratio); | |
1406 | cl2->SetType(10); | |
1407 | pairs[ip][j]=10; | |
1408 | if ((pos[ip].GetNd()+neg[j].GetNd())>6){ //multi cluster | |
1409 | cl2->SetType(11); | |
1410 | pairs[ip][j]=11; | |
1411 | } | |
1412 | cused1[ip]++; | |
1413 | cused2[j]++; | |
1414 | } | |
1415 | ||
1416 | } | |
1417 | ||
1418 | // | |
c391f9d9 | 1419 | for (Int_t i=0; i<np; i++) { |
c391f9d9 | 1420 | Float_t ybest=1000,zbest=1000,qbest=0; |
1421 | Float_t yp=pos[i].GetY()*fYpitchSSD; | |
e43c066c | 1422 | if (pos[i].GetQ()<3) continue; |
c391f9d9 | 1423 | for (Int_t j=0; j<nn; j++) { |
e43c066c | 1424 | // for (Int_t di = 0;di<cpositive[i];di++){ |
1425 | // Int_t j = negativepair[10*i+di]; | |
1426 | if (neg[j].GetQ()<3) continue; | |
1427 | if (cused2[j]||cused1[i]) continue; | |
1428 | if (pairs[i][j]>0 &&pairs[i][j]<100) continue; | |
1429 | ratio = (pos[i].GetQ()-neg[j].GetQ())/(pos[i].GetQ()+neg[j].GetQ()); | |
c391f9d9 | 1430 | Float_t yn=neg[j].GetY()*fYpitchSSD; |
1431 | Float_t zt=(2*fHlSSD*tanp + yp - yn)/(tann+tanp); | |
1432 | Float_t yt=yn + tann*zt; | |
1433 | zt-=fHlSSD; yt-=fHwSSD; | |
1434 | if (TMath::Abs(yt)<fHwSSD+0.01) | |
e43c066c | 1435 | if (TMath::Abs(zt)<fHlSSD+0.01*(neg[j].GetNd()+pos[i].GetNd())) { |
c391f9d9 | 1436 | ybest=yt; zbest=zt; |
1437 | qbest=0.5*(pos[i].GetQ()+neg[j].GetQ()); | |
f363d377 | 1438 | lp[0]=-(-ybest+fYshift[fI]); |
1439 | lp[1]= -zbest+fZshift[fI]; | |
e43c066c | 1440 | lp[2]=0.0025*0.0025; //SigmaY2 |
1441 | lp[3]=0.110*0.110; //SigmaZ2 | |
1442 | ||
c391f9d9 | 1443 | lp[4]=qbest; //Q |
babd135a | 1444 | for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2; |
e43c066c | 1445 | for (Int_t ilab=0;ilab<3;ilab++){ |
1446 | milab[ilab] = pos[i].GetLabel(ilab); | |
1447 | milab[ilab+3] = neg[j].GetLabel(ilab); | |
1448 | } | |
babd135a | 1449 | // |
babd135a | 1450 | CheckLabels2(milab); |
1451 | milab[3]=(((i<<10) + j)<<10) + idet; // pos|neg|det | |
e43c066c | 1452 | Int_t info[3] = {pos[i].GetNd(),neg[j].GetNd(),fNlayer[fI]}; |
1453 | AliITSclusterV2 * cl2 = new (cl[ncl]) AliITSclusterV2(milab,lp,info); | |
1454 | ncl++; | |
1455 | cl2->SetChargeRatio(ratio); | |
1456 | cl2->SetType(100+cpositive[j]+cnegative[i]); | |
1457 | //cl2->SetType(0); | |
1458 | /* | |
1459 | if (pairs[i][j]<100){ | |
1460 | printf("problem:- %d\n", pairs[i][j]); | |
1461 | } | |
1462 | if (cnegative[i]<2&&cpositive[j]<2){ | |
1463 | printf("problem:- %d\n", pairs[i][j]); | |
1464 | } | |
1465 | */ | |
c391f9d9 | 1466 | } |
1467 | } | |
c391f9d9 | 1468 | } |
a06cb96e | 1469 | |
0774985e | 1470 | // for (Int_t i=0; i<1000; i++) delete [] pairs[i]; |
1471 | // delete [] pairs; | |
a06cb96e | 1472 | |
c391f9d9 | 1473 | } |
1474 | ||
e43c066c | 1475 | |
7f4044f0 | 1476 | void AliITSclustererV2:: |
e43c066c | 1477 | FindClustersSSD(const TClonesArray *alldigits, TClonesArray *clusters) { |
7f4044f0 | 1478 | //------------------------------------------------------------ |
1479 | // Actual SSD cluster finder | |
1480 | //------------------------------------------------------------ | |
e43c066c | 1481 | Int_t smaxall=alldigits->GetEntriesFast(); |
1482 | if (smaxall==0) return; | |
1483 | TObjArray *digits = new TObjArray; | |
1484 | for (Int_t i=0;i<smaxall; i++){ | |
1485 | AliITSdigitSSD *d=(AliITSdigitSSD*)alldigits->UncheckedAt(i); | |
1486 | if (d->GetSignal()<3) continue; | |
1487 | digits->AddLast(d); | |
1488 | } | |
1489 | Int_t smax = digits->GetEntriesFast(); | |
7f4044f0 | 1490 | if (smax==0) return; |
e43c066c | 1491 | |
7f4044f0 | 1492 | const Int_t MAX=1000; |
1493 | Int_t np=0, nn=0; | |
1494 | Ali1Dcluster pos[MAX], neg[MAX]; | |
1495 | Float_t y=0., q=0., qmax=0.; | |
1496 | Int_t lab[4]={-2,-2,-2,-2}; | |
e43c066c | 1497 | |
7f4044f0 | 1498 | AliITSdigitSSD *d=(AliITSdigitSSD*)digits->UncheckedAt(0); |
1499 | q += d->GetSignal(); | |
1500 | y += d->GetCoord2()*d->GetSignal(); | |
1501 | qmax=d->GetSignal(); | |
1502 | lab[0]=d->GetTrack(0); lab[1]=d->GetTrack(1); lab[2]=d->GetTrack(2); | |
1503 | Int_t curr=d->GetCoord2(); | |
1504 | Int_t flag=d->GetCoord1(); | |
1505 | Int_t *n=&nn; | |
1506 | Ali1Dcluster *c=neg; | |
c391f9d9 | 1507 | Int_t nd=1; |
e43c066c | 1508 | Int_t milab[10]; |
1509 | for (Int_t ilab=0;ilab<10;ilab++){ | |
1510 | milab[ilab]=-2; | |
1511 | } | |
1512 | milab[0]=d->GetTrack(0); milab[1]=d->GetTrack(1); milab[2]=d->GetTrack(2); | |
1513 | ||
7f4044f0 | 1514 | for (Int_t s=1; s<smax; s++) { |
e43c066c | 1515 | d=(AliITSdigitSSD*)digits->UncheckedAt(s); |
7f4044f0 | 1516 | Int_t strip=d->GetCoord2(); |
1517 | if ((strip-curr) > 1 || flag!=d->GetCoord1()) { | |
1518 | c[*n].SetY(y/q); | |
1519 | c[*n].SetQ(q); | |
1520 | c[*n].SetNd(nd); | |
e43c066c | 1521 | CheckLabels2(milab); |
1522 | c[*n].SetLabels(milab); | |
7f4044f0 | 1523 | //Split suspiciously big cluster |
e43c066c | 1524 | /* |
1525 | if (nd>10&&nd<16){ | |
1526 | c[*n].SetY(y/q-0.3*nd); | |
1527 | c[*n].SetQ(0.5*q); | |
1528 | (*n)++; | |
1529 | if (*n==MAX) { | |
1530 | Error("FindClustersSSD","Too many 1D clusters !"); | |
7f4044f0 | 1531 | return; |
e43c066c | 1532 | } |
1533 | c[*n].SetY(y/q-0.0*nd); | |
1534 | c[*n].SetQ(0.5*q); | |
1535 | c[*n].SetNd(nd); | |
1536 | (*n)++; | |
1537 | if (*n==MAX) { | |
1538 | Error("FindClustersSSD","Too many 1D clusters !"); | |
1539 | return; | |
1540 | } | |
1541 | // | |
1542 | c[*n].SetY(y/q+0.3*nd); | |
1543 | c[*n].SetQ(0.5*q); | |
1544 | c[*n].SetNd(nd); | |
1545 | c[*n].SetLabels(milab); | |
1546 | } | |
1547 | else{ | |
1548 | */ | |
1549 | if (nd>4&&nd<25) { | |
1550 | c[*n].SetY(y/q-0.25*nd); | |
1551 | c[*n].SetQ(0.5*q); | |
1552 | (*n)++; | |
1553 | if (*n==MAX) { | |
1554 | Error("FindClustersSSD","Too many 1D clusters !"); | |
1555 | return; | |
1556 | } | |
1557 | c[*n].SetY(y/q+0.25*nd); | |
1558 | c[*n].SetQ(0.5*q); | |
1559 | c[*n].SetNd(nd); | |
1560 | c[*n].SetLabels(milab); | |
1561 | } | |
7f4044f0 | 1562 | (*n)++; |
1563 | if (*n==MAX) { | |
c630aafd | 1564 | Error("FindClustersSSD","Too many 1D clusters !"); |
7f4044f0 | 1565 | return; |
1566 | } | |
1567 | y=q=qmax=0.; | |
1568 | nd=0; | |
1569 | lab[0]=lab[1]=lab[2]=-2; | |
e43c066c | 1570 | // |
1571 | for (Int_t ilab=0;ilab<10;ilab++){ | |
1572 | milab[ilab]=-2; | |
1573 | } | |
1574 | // | |
7f4044f0 | 1575 | if (flag!=d->GetCoord1()) { n=&np; c=pos; } |
1576 | } | |
1577 | flag=d->GetCoord1(); | |
1578 | q += d->GetSignal(); | |
1579 | y += d->GetCoord2()*d->GetSignal(); | |
1580 | nd++; | |
1581 | if (d->GetSignal()>qmax) { | |
1582 | qmax=d->GetSignal(); | |
1583 | lab[0]=d->GetTrack(0); lab[1]=d->GetTrack(1); lab[2]=d->GetTrack(2); | |
1584 | } | |
e43c066c | 1585 | for (Int_t ilab=0;ilab<10;ilab++) { |
1586 | if (d->GetTrack(ilab)>=0) AddLabel(milab, (d->GetTrack(ilab))); | |
1587 | } | |
7f4044f0 | 1588 | curr=strip; |
1589 | } | |
1590 | c[*n].SetY(y/q); | |
1591 | c[*n].SetQ(q); | |
1592 | c[*n].SetNd(nd); | |
1593 | c[*n].SetLabels(lab); | |
1594 | //Split suspiciously big cluster | |
0e986ef9 | 1595 | if (nd>4 && nd<25) { |
1596 | c[*n].SetY(y/q-0.25*nd); | |
7f4044f0 | 1597 | c[*n].SetQ(0.5*q); |
1598 | (*n)++; | |
1599 | if (*n==MAX) { | |
c630aafd | 1600 | Error("FindClustersSSD","Too many 1D clusters !"); |
7f4044f0 | 1601 | return; |
1602 | } | |
0e986ef9 | 1603 | c[*n].SetY(y/q+0.25*nd); |
7f4044f0 | 1604 | c[*n].SetQ(0.5*q); |
1605 | c[*n].SetNd(nd); | |
1606 | c[*n].SetLabels(lab); | |
1607 | } | |
1608 | (*n)++; | |
1609 | if (*n==MAX) { | |
c630aafd | 1610 | Error("FindClustersSSD","Too many 1D clusters !"); |
7f4044f0 | 1611 | return; |
1612 | } | |
1613 | ||
c391f9d9 | 1614 | FindClustersSSD(neg, nn, pos, np, clusters); |
1615 | } | |
7f4044f0 | 1616 | |
c391f9d9 | 1617 | void AliITSclustererV2::FindClustersSSD(AliITSRawStream* input, |
1618 | TClonesArray** clusters) | |
1619 | { | |
1620 | //------------------------------------------------------------ | |
1621 | // Actual SSD cluster finder for raw data | |
1622 | //------------------------------------------------------------ | |
1623 | Int_t nClustersSSD = 0; | |
1624 | const Int_t MAX = 1000; | |
1625 | Ali1Dcluster clusters1D[2][MAX]; | |
1626 | Int_t nClusters[2] = {0, 0}; | |
0e986ef9 | 1627 | Int_t lab[3]={-2,-2,-2}; |
c391f9d9 | 1628 | Float_t q = 0.; |
1629 | Float_t y = 0.; | |
1630 | Int_t nDigits = 0; | |
1631 | Int_t prevStrip = -1; | |
1632 | Int_t prevFlag = -1; | |
0e986ef9 | 1633 | Int_t prevModule = -1; |
c391f9d9 | 1634 | |
1635 | // read raw data input stream | |
1636 | while (kTRUE) { | |
1637 | Bool_t next = input->Next(); | |
1638 | ||
0e986ef9 | 1639 | if(input->GetSignal()<3 && next) continue; |
c391f9d9 | 1640 | // check if a new cluster starts |
1641 | Int_t strip = input->GetCoord2(); | |
1642 | Int_t flag = input->GetCoord1(); | |
0e986ef9 | 1643 | if ((!next || (input->GetModuleID() != prevModule)|| |
c391f9d9 | 1644 | (strip-prevStrip > 1) || (flag != prevFlag)) && |
1645 | (nDigits > 0)) { | |
1646 | if (nClusters[prevFlag] == MAX) { | |
1647 | Error("FindClustersSSD", "Too many 1D clusters !"); | |
1648 | return; | |
1649 | } | |
1650 | Ali1Dcluster& cluster = clusters1D[prevFlag][nClusters[prevFlag]++]; | |
1651 | cluster.SetY(y/q); | |
1652 | cluster.SetQ(q); | |
1653 | cluster.SetNd(nDigits); | |
0e986ef9 | 1654 | cluster.SetLabels(lab); |
c391f9d9 | 1655 | |
1656 | //Split suspiciously big cluster | |
0e986ef9 | 1657 | if (nDigits > 4&&nDigits < 25) { |
1658 | cluster.SetY(y/q - 0.25*nDigits); | |
c391f9d9 | 1659 | cluster.SetQ(0.5*q); |
1660 | if (nClusters[prevFlag] == MAX) { | |
1661 | Error("FindClustersSSD", "Too many 1D clusters !"); | |
1662 | return; | |
1663 | } | |
1664 | Ali1Dcluster& cluster2 = clusters1D[prevFlag][nClusters[prevFlag]++]; | |
0e986ef9 | 1665 | cluster2.SetY(y/q + 0.25*nDigits); |
c391f9d9 | 1666 | cluster2.SetQ(0.5*q); |
1667 | cluster2.SetNd(nDigits); | |
0e986ef9 | 1668 | cluster2.SetLabels(lab); |
7f4044f0 | 1669 | } |
c391f9d9 | 1670 | y = q = 0.; |
1671 | nDigits = 0; | |
7f4044f0 | 1672 | } |
c391f9d9 | 1673 | |
0e986ef9 | 1674 | if (!next || (input->GetModuleID() != prevModule)) { |
1675 | Int_t iModule = prevModule; | |
c391f9d9 | 1676 | |
1677 | // when all data from a module was read, search for clusters | |
1678 | if (prevFlag >= 0) { | |
1679 | clusters[iModule] = new TClonesArray("AliITSclusterV2"); | |
1680 | fI = iModule; | |
1681 | FindClustersSSD(&clusters1D[0][0], nClusters[0], | |
1682 | &clusters1D[1][0], nClusters[1], clusters[iModule]); | |
1683 | Int_t nClusters = clusters[iModule]->GetEntriesFast(); | |
1684 | nClustersSSD += nClusters; | |
1685 | } | |
1686 | ||
1687 | if (!next) break; | |
1688 | nClusters[0] = nClusters[1] = 0; | |
1689 | y = q = 0.; | |
1690 | nDigits = 0; | |
7f4044f0 | 1691 | } |
c391f9d9 | 1692 | |
1693 | // add digit to current cluster | |
1694 | q += input->GetSignal(); | |
1695 | y += strip * input->GetSignal(); | |
1696 | nDigits++; | |
1697 | prevStrip = strip; | |
1698 | prevFlag = flag; | |
0e986ef9 | 1699 | prevModule = input->GetModuleID(); |
1700 | ||
7f4044f0 | 1701 | } |
c391f9d9 | 1702 | |
1703 | Info("FindClustersSSD", "found clusters in ITS SSD: %d", nClustersSSD); | |
7f4044f0 | 1704 | } |
1705 | ||
1706 | #else //V1 | |
1707 | ||
1708 | #include "AliITSDetType.h" | |
1709 | ||
1710 | #include "AliITSsegmentationSPD.h" | |
1711 | #include "AliITSClusterFinderSPD.h" | |
1712 | ||
1713 | #include "AliITSresponseSDD.h" | |
1714 | #include "AliITSsegmentationSDD.h" | |
1715 | #include "AliITSClusterFinderSDD.h" | |
1716 | ||
1717 | #include "AliITSsegmentationSSD.h" | |
1718 | #include "AliITSClusterFinderSSD.h" | |
1719 | ||
1720 | ||
1721 | void AliITSclustererV2:: | |
1722 | FindClustersSPD(const TClonesArray *digits, TClonesArray *clusters) { | |
1723 | //------------------------------------------------------------ | |
1724 | // Actual SPD cluster finding based on AliITSClusterFinderSPD | |
1725 | //------------------------------------------------------------ | |
1726 | static AliITS *its=(AliITS*)gAlice->GetModule("ITS"); | |
1727 | static TClonesArray *points=its->RecPoints(); | |
1728 | static AliITSsegmentationSPD *seg= | |
1729 | (AliITSsegmentationSPD *)its->DetType(0)->GetSegmentationModel(); | |
1730 | static AliITSClusterFinderSPD cf(seg, (TClonesArray*)digits, points); | |
1731 | ||
1732 | cf.FindRawClusters(fI); | |
1cca57bf | 1733 | RecPoints2Clusters(points, fI, clusters); |
7f4044f0 | 1734 | its->ResetRecPoints(); |
1735 | ||
1736 | } | |
1737 | ||
1738 | void AliITSclustererV2:: | |
1739 | FindClustersSDD(const TClonesArray *digits, TClonesArray *clusters) { | |
1740 | //------------------------------------------------------------ | |
1741 | // Actual SDD cluster finding based on AliITSClusterFinderSDD | |
1742 | //------------------------------------------------------------ | |
1743 | static AliITS *its=(AliITS*)gAlice->GetModule("ITS"); | |
1744 | static TClonesArray *points=its->RecPoints(); | |
1745 | static AliITSresponseSDD *resp= | |
1746 | (AliITSresponseSDD *)its->DetType(1)->GetResponseModel(); | |
1747 | static AliITSsegmentationSDD *seg= | |
1748 | (AliITSsegmentationSDD *)its->DetType(1)->GetSegmentationModel(); | |
1749 | static AliITSClusterFinderSDD | |
1750 | cf(seg,resp,(TClonesArray*)digits,its->ClustersAddress(1)); | |
1751 | ||
1752 | cf.FindRawClusters(fI); | |
1cca57bf | 1753 | Int_t nc=points->GetEntriesFast(); |
1754 | while (nc--) { //To be consistent with the SSD cluster charges | |
1755 | AliITSRecPoint *p=(AliITSRecPoint*)points->UncheckedAt(nc); | |
e3d91d99 | 1756 | p->SetQ(p->GetQ()/12.); |
1cca57bf | 1757 | } |
1758 | RecPoints2Clusters(points, fI, clusters); | |
7f4044f0 | 1759 | its->ResetClusters(1); |
1760 | its->ResetRecPoints(); | |
1761 | ||
1762 | } | |
1763 | ||
1764 | void AliITSclustererV2:: | |
1765 | FindClustersSSD(const TClonesArray *digits, TClonesArray *clusters) { | |
1766 | //------------------------------------------------------------ | |
1767 | // Actual SSD cluster finding based on AliITSClusterFinderSSD | |
1768 | //------------------------------------------------------------ | |
1769 | static AliITS *its=(AliITS*)gAlice->GetModule("ITS"); | |
1770 | static TClonesArray *points=its->RecPoints(); | |
1771 | static AliITSsegmentationSSD *seg= | |
1772 | (AliITSsegmentationSSD *)its->DetType(2)->GetSegmentationModel(); | |
1773 | static AliITSClusterFinderSSD cf(seg,(TClonesArray*)digits); | |
1774 | ||
1775 | cf.FindRawClusters(fI); | |
1cca57bf | 1776 | RecPoints2Clusters(points, fI, clusters); |
7f4044f0 | 1777 | its->ResetRecPoints(); |
1778 | ||
1779 | } | |
1780 | ||
1781 | #endif |