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