]> git.uio.no Git - u/mrichter/AliRoot.git/blame - ITS/AliITSClusterFinderV2SSD.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
changes to be compliant with Eff C++ rules
[u/mrichter/AliRoot.git] / ITS / AliITSClusterFinderV2SSD.cxx
CommitLineData
04366a57 1/**************************************************************************
2 * Copyright(c) 1998-2003, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
15////////////////////////////////////////////////////////////////////////////
16// Implementation of the ITS clusterer V2 class //
17// //
18// Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch //
19// //
20///////////////////////////////////////////////////////////////////////////
21
04366a57 22
23#include "AliITSClusterFinderV2SSD.h"
00a7cc50 24#include "AliITSRecPoint.h"
7d62fb64 25#include "AliITSDetTypeRec.h"
04366a57 26#include "AliRawReader.h"
27#include "AliITSRawStreamSSD.h"
04366a57 28#include <TClonesArray.h>
04366a57 29#include "AliITSdigitSSD.h"
30
31ClassImp(AliITSClusterFinderV2SSD)
32
33
e56160b8 34AliITSClusterFinderV2SSD::AliITSClusterFinderV2SSD(AliITSDetTypeRec* dettyp):AliITSClusterFinderV2(dettyp),
35fLastSSD1(0),
36fYpitchSSD(0.0095),
37fHwSSD(3.65),
38fHlSSD(2.00),
39fTanP(0.0275),
40fTanN(0.0075){
04366a57 41
42 //Default constructor
43
8ba39da9 44 fLastSSD1=fDetTypeRec->GetITSgeom()->GetModuleIndex(6,1,1)-1;
04366a57 45
46}
47
48
49void AliITSClusterFinderV2SSD::FindRawClusters(Int_t mod){
50
51 //Find clusters V2
52 SetModule(mod);
53 FindClustersSSD(fDigits);
54
55}
56
57void AliITSClusterFinderV2SSD::FindClustersSSD(TClonesArray *alldigits) {
58 //------------------------------------------------------------
59 // Actual SSD cluster finder
60 //------------------------------------------------------------
61 Int_t smaxall=alldigits->GetEntriesFast();
62 if (smaxall==0) return;
63 TObjArray *digits = new TObjArray;
64 for (Int_t i=0;i<smaxall; i++){
65 AliITSdigitSSD *d=(AliITSdigitSSD*)alldigits->UncheckedAt(i);
66 if (d->GetSignal()<3) continue;
67 digits->AddLast(d);
68 }
69 Int_t smax = digits->GetEntriesFast();
70 if (smax==0) return;
71
72 const Int_t kMax=1000;
73 Int_t np=0, nn=0;
74 Ali1Dcluster pos[kMax], neg[kMax];
75 Float_t y=0., q=0., qmax=0.;
76 Int_t lab[4]={-2,-2,-2,-2};
77
78 AliITSdigitSSD *d=(AliITSdigitSSD*)digits->UncheckedAt(0);
79 q += d->GetSignal();
80 y += d->GetCoord2()*d->GetSignal();
81 qmax=d->GetSignal();
82 lab[0]=d->GetTrack(0); lab[1]=d->GetTrack(1); lab[2]=d->GetTrack(2);
83 Int_t curr=d->GetCoord2();
84 Int_t flag=d->GetCoord1();
85 Int_t *n=&nn;
86 Ali1Dcluster *c=neg;
87 Int_t nd=1;
88 Int_t milab[10];
89 for (Int_t ilab=0;ilab<10;ilab++){
90 milab[ilab]=-2;
91 }
92 milab[0]=d->GetTrack(0); milab[1]=d->GetTrack(1); milab[2]=d->GetTrack(2);
93
94 for (Int_t s=1; s<smax; s++) {
95 d=(AliITSdigitSSD*)digits->UncheckedAt(s);
96 Int_t strip=d->GetCoord2();
97 if ((strip-curr) > 1 || flag!=d->GetCoord1()) {
98 c[*n].SetY(y/q);
99 c[*n].SetQ(q);
100 c[*n].SetNd(nd);
101 CheckLabels2(milab);
102 c[*n].SetLabels(milab);
103 //Split suspiciously big cluster
104 /*
105 if (nd>10&&nd<16){
106 c[*n].SetY(y/q-0.3*nd);
107 c[*n].SetQ(0.5*q);
108 (*n)++;
109 if (*n==MAX) {
110 Error("FindClustersSSD","Too many 1D clusters !");
111 return;
112 }
113 c[*n].SetY(y/q-0.0*nd);
114 c[*n].SetQ(0.5*q);
115 c[*n].SetNd(nd);
116 (*n)++;
117 if (*n==MAX) {
118 Error("FindClustersSSD","Too many 1D clusters !");
119 return;
120 }
121 //
122 c[*n].SetY(y/q+0.3*nd);
123 c[*n].SetQ(0.5*q);
124 c[*n].SetNd(nd);
125 c[*n].SetLabels(milab);
126 }
127 else{
128 */
129 if (nd>4&&nd<25) {
130 c[*n].SetY(y/q-0.25*nd);
131 c[*n].SetQ(0.5*q);
132 (*n)++;
133 if (*n==kMax) {
134 Error("FindClustersSSD","Too many 1D clusters !");
135 return;
136 }
137 c[*n].SetY(y/q+0.25*nd);
138 c[*n].SetQ(0.5*q);
139 c[*n].SetNd(nd);
140 c[*n].SetLabels(milab);
141 }
142 (*n)++;
143 if (*n==kMax) {
144 Error("FindClustersSSD","Too many 1D clusters !");
145 return;
146 }
147 y=q=qmax=0.;
148 nd=0;
149 lab[0]=lab[1]=lab[2]=-2;
150 //
151 for (Int_t ilab=0;ilab<10;ilab++){
152 milab[ilab]=-2;
153 }
154 //
155 if (flag!=d->GetCoord1()) { n=&np; c=pos; }
156 }
157 flag=d->GetCoord1();
158 q += d->GetSignal();
159 y += d->GetCoord2()*d->GetSignal();
160 nd++;
161 if (d->GetSignal()>qmax) {
162 qmax=d->GetSignal();
163 lab[0]=d->GetTrack(0); lab[1]=d->GetTrack(1); lab[2]=d->GetTrack(2);
164 }
165 for (Int_t ilab=0;ilab<10;ilab++) {
166 if (d->GetTrack(ilab)>=0) AddLabel(milab, (d->GetTrack(ilab)));
167 }
168 curr=strip;
169 }
170 c[*n].SetY(y/q);
171 c[*n].SetQ(q);
172 c[*n].SetNd(nd);
173 c[*n].SetLabels(lab);
174 //Split suspiciously big cluster
175 if (nd>4 && nd<25) {
176 c[*n].SetY(y/q-0.25*nd);
177 c[*n].SetQ(0.5*q);
178 (*n)++;
179 if (*n==kMax) {
180 Error("FindClustersSSD","Too many 1D clusters !");
181 return;
182 }
183 c[*n].SetY(y/q+0.25*nd);
184 c[*n].SetQ(0.5*q);
185 c[*n].SetNd(nd);
186 c[*n].SetLabels(lab);
187 }
188 (*n)++;
189 if (*n==kMax) {
190 Error("FindClustersSSD","Too many 1D clusters !");
191 return;
192 }
193
194 FindClustersSSD(neg, nn, pos, np);
195}
196
197
198void AliITSClusterFinderV2SSD::RawdataToClusters(AliRawReader* rawReader,TClonesArray** clusters){
199
200 //------------------------------------------------------------
201 // This function creates ITS clusters from raw data
202 //------------------------------------------------------------
203 rawReader->Reset();
204 AliITSRawStreamSSD inputSSD(rawReader);
205 FindClustersSSD(&inputSSD,clusters);
206
207}
208
209void AliITSClusterFinderV2SSD::FindClustersSSD(AliITSRawStream* input,
210 TClonesArray** clusters)
211{
212 //------------------------------------------------------------
213 // Actual SSD cluster finder for raw data
214 //------------------------------------------------------------
215 Int_t nClustersSSD = 0;
216 const Int_t kMax = 1000;
217 Ali1Dcluster clusters1D[2][kMax];
218 Int_t nClusters[2] = {0, 0};
219 Int_t lab[3]={-2,-2,-2};
220 Float_t q = 0.;
221 Float_t y = 0.;
222 Int_t nDigits = 0;
223 Int_t prevStrip = -1;
224 Int_t prevFlag = -1;
225 Int_t prevModule = -1;
226
227 // read raw data input stream
228 while (kTRUE) {
229 Bool_t next = input->Next();
230
231 if(input->GetSignal()<3 && next) continue;
232 // check if a new cluster starts
233 Int_t strip = input->GetCoord2();
234 Int_t flag = input->GetCoord1();
235 if ((!next || (input->GetModuleID() != prevModule)||
236 (strip-prevStrip > 1) || (flag != prevFlag)) &&
237 (nDigits > 0)) {
238 if (nClusters[prevFlag] == kMax) {
239 Error("FindClustersSSD", "Too many 1D clusters !");
240 return;
241 }
242 Ali1Dcluster& cluster = clusters1D[prevFlag][nClusters[prevFlag]++];
243 cluster.SetY(y/q);
244 cluster.SetQ(q);
245 cluster.SetNd(nDigits);
246 cluster.SetLabels(lab);
247
248 //Split suspiciously big cluster
249 if (nDigits > 4&&nDigits < 25) {
250 cluster.SetY(y/q - 0.25*nDigits);
251 cluster.SetQ(0.5*q);
252 if (nClusters[prevFlag] == kMax) {
253 Error("FindClustersSSD", "Too many 1D clusters !");
254 return;
255 }
256 Ali1Dcluster& cluster2 = clusters1D[prevFlag][nClusters[prevFlag]++];
257 cluster2.SetY(y/q + 0.25*nDigits);
258 cluster2.SetQ(0.5*q);
259 cluster2.SetNd(nDigits);
260 cluster2.SetLabels(lab);
261 }
262 y = q = 0.;
263 nDigits = 0;
264 }
265
266 if (!next || (input->GetModuleID() != prevModule)) {
267 Int_t iModule = prevModule;
268
269 // when all data from a module was read, search for clusters
270 if (prevFlag >= 0) {
00a7cc50 271 clusters[iModule] = new TClonesArray("AliITSRecPoint");
04366a57 272 fModule = iModule;
273 FindClustersSSD(&clusters1D[0][0], nClusters[0],
274 &clusters1D[1][0], nClusters[1], clusters[iModule]);
275 Int_t nClusters = clusters[iModule]->GetEntriesFast();
276 nClustersSSD += nClusters;
277 }
278
279 if (!next) break;
280 nClusters[0] = nClusters[1] = 0;
281 y = q = 0.;
282 nDigits = 0;
283 }
284
285 // add digit to current cluster
286 q += input->GetSignal();
287 y += strip * input->GetSignal();
288 nDigits++;
289 prevStrip = strip;
290 prevFlag = flag;
291 prevModule = input->GetModuleID();
292
293 }
294
295 Info("FindClustersSSD", "found clusters in ITS SSD: %d", nClustersSSD);
296}
297
298void AliITSClusterFinderV2SSD::
299FindClustersSSD(Ali1Dcluster* neg, Int_t nn,
300 Ali1Dcluster* pos, Int_t np,
301 TClonesArray *clusters) {
302 //------------------------------------------------------------
303 // Actual SSD cluster finder
304 //------------------------------------------------------------
305 TClonesArray &cl=*clusters;
306 //
307 Float_t tanp=fTanP, tann=fTanN;
308 if (fModule>fLastSSD1) {tann=fTanP; tanp=fTanN;}
309 Int_t idet=fNdet[fModule];
310 Int_t ncl=0;
311 //
312 Int_t negativepair[30000];
313 Int_t cnegative[3000];
314 Int_t cused1[3000];
315 Int_t positivepair[30000];
316 Int_t cpositive[3000];
317 Int_t cused2[3000];
318 for (Int_t i=0;i<3000;i++) {cnegative[i]=0; cused1[i]=0;}
319 for (Int_t i=0;i<3000;i++) {cpositive[i]=0; cused2[i]=0;}
78538bfe 320 for (Int_t i=0;i<30000;i++) {negativepair[i]=0; positivepair[i]=0;}
04366a57 321 static Short_t pairs[1000][1000];
322 memset(pairs,0,sizeof(Short_t)*1000000);
323// Short_t ** pairs = new Short_t*[1000];
324// for (Int_t i=0; i<1000; i++) {
325// pairs[i] = new Short_t[1000];
326// memset(pairs[i],0,sizeof(Short_t)*1000);
327// }
328 //
329 // find available pairs
330 //
331 for (Int_t i=0; i<np; i++) {
332 Float_t yp=pos[i].GetY()*fYpitchSSD;
333 if (pos[i].GetQ()<3) continue;
334 for (Int_t j=0; j<nn; j++) {
335 if (neg[j].GetQ()<3) continue;
336 Float_t yn=neg[j].GetY()*fYpitchSSD;
337 Float_t zt=(2*fHlSSD*tanp + yp - yn)/(tann+tanp);
338 Float_t yt=yn + tann*zt;
339 zt-=fHlSSD; yt-=fHwSSD;
340 if (TMath::Abs(yt)<fHwSSD+0.01)
341 if (TMath::Abs(zt)<fHlSSD+0.01*(neg[j].GetNd()+pos[i].GetNd())) {
342 negativepair[i*10+cnegative[i]] =j; //index
343 positivepair[j*10+cpositive[j]] =i;
344 cnegative[i]++; //counters
345 cpositive[j]++;
346 pairs[i][j]=100;
347 }
348 }
349 }
350 //
351 for (Int_t i=0; i<np; i++) {
352 Float_t yp=pos[i].GetY()*fYpitchSSD;
353 if (pos[i].GetQ()<3) continue;
354 for (Int_t j=0; j<nn; j++) {
355 if (neg[j].GetQ()<3) continue;
356 if (cpositive[j]&&cnegative[i]) continue;
357 Float_t yn=neg[j].GetY()*fYpitchSSD;
358 Float_t zt=(2*fHlSSD*tanp + yp - yn)/(tann+tanp);
359 Float_t yt=yn + tann*zt;
360 zt-=fHlSSD; yt-=fHwSSD;
361 if (TMath::Abs(yt)<fHwSSD+0.1)
362 if (TMath::Abs(zt)<fHlSSD+0.15) {
363 if (cnegative[i]==0) pos[i].SetNd(100); // not available pair
364 if (cpositive[j]==0) neg[j].SetNd(100); // not available pair
365 negativepair[i*10+cnegative[i]] =j; //index
366 positivepair[j*10+cpositive[j]] =i;
367 cnegative[i]++; //counters
368 cpositive[j]++;
369 pairs[i][j]=100;
370 }
371 }
372 }
373 //
374 Float_t lp[5];
375 Int_t milab[10];
376 Double_t ratio;
377
378 //
379 // sign gold tracks
380 //
381 for (Int_t ip=0;ip<np;ip++){
382 Float_t ybest=1000,zbest=1000,qbest=0;
383 //
384 // select gold clusters
385 if ( (cnegative[ip]==1) && cpositive[negativepair[10*ip]]==1){
386 Float_t yp=pos[ip].GetY()*fYpitchSSD;
387 Int_t j = negativepair[10*ip];
388 ratio = (pos[ip].GetQ()-neg[j].GetQ())/(pos[ip].GetQ()+neg[j].GetQ());
389 //
390 Float_t yn=neg[j].GetY()*fYpitchSSD;
391 Float_t zt=(2*fHlSSD*tanp + yp - yn)/(tann+tanp);
392 Float_t yt=yn + tann*zt;
393 zt-=fHlSSD; yt-=fHwSSD;
394 ybest=yt; zbest=zt;
395 qbest=0.5*(pos[ip].GetQ()+neg[j].GetQ());
396 lp[0]=-(-ybest+fYshift[fModule]);
397 lp[1]= -zbest+fZshift[fModule];
398 lp[2]=0.0025*0.0025; //SigmaY2
399 lp[3]=0.110*0.110; //SigmaZ2
400
401 lp[4]=qbest; //Q
402 for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2;
403 for (Int_t ilab=0;ilab<3;ilab++){
404 milab[ilab] = pos[ip].GetLabel(ilab);
405 milab[ilab+3] = neg[j].GetLabel(ilab);
406 }
407 //
408 CheckLabels2(milab);
409 milab[3]=(((ip<<10) + j)<<10) + idet; // pos|neg|det
410 Int_t info[3] = {pos[ip].GetNd(),neg[j].GetNd(),fNlayer[fModule]};
00a7cc50 411 AliITSRecPoint * cl2;
412 if(clusters){
413 cl2 = new (cl[ncl]) AliITSRecPoint(fModule,fDetTypeRec->GetITSgeom(),milab,lp,info);
414 cl2->SetChargeRatio(ratio);
415 cl2->SetType(1);
416 pairs[ip][j]=1;
417 if ((pos[ip].GetNd()+neg[j].GetNd())>6){ //multi cluster
418 cl2->SetType(2);
419 pairs[ip][j]=2;
420 }
421 cused1[ip]++;
422 cused2[j]++;
423
424 }
04366a57 425 else{
00a7cc50 426 cl2 = new AliITSRecPoint(fModule,fDetTypeRec->GetITSgeom(),milab,lp,info);
427 cl2->SetChargeRatio(ratio);
428 cl2->SetType(1);
429 pairs[ip][j]=1;
430 if ((pos[ip].GetNd()+neg[j].GetNd())>6){ //multi cluster
431 cl2->SetType(2);
432 pairs[ip][j]=2;
433 }
434 cused1[ip]++;
435 cused2[j]++;
436 fDetTypeRec->AddRecPoint(*cl2);
04366a57 437 }
438 ncl++;
04366a57 439 }
440 }
441
442 for (Int_t ip=0;ip<np;ip++){
443 Float_t ybest=1000,zbest=1000,qbest=0;
444 //
445 //
446 // select "silber" cluster
447 if ( cnegative[ip]==1 && cpositive[negativepair[10*ip]]==2){
448 Int_t in = negativepair[10*ip];
449 Int_t ip2 = positivepair[10*in];
450 if (ip2==ip) ip2 = positivepair[10*in+1];
451 Float_t pcharge = pos[ip].GetQ()+pos[ip2].GetQ();
452 if (TMath::Abs(pcharge-neg[in].GetQ())<10){
453 //
454 // add first pair
455 if (pairs[ip][in]==100){ //
456 Float_t yp=pos[ip].GetY()*fYpitchSSD;
457 Float_t yn=neg[in].GetY()*fYpitchSSD;
458 Float_t zt=(2*fHlSSD*tanp + yp - yn)/(tann+tanp);
459 Float_t yt=yn + tann*zt;
460 zt-=fHlSSD; yt-=fHwSSD;
461 ybest =yt; zbest=zt;
462 qbest =pos[ip].GetQ();
463 lp[0]=-(-ybest+fYshift[fModule]);
464 lp[1]= -zbest+fZshift[fModule];
465 lp[2]=0.0025*0.0025; //SigmaY2
466 lp[3]=0.110*0.110; //SigmaZ2
467
468 lp[4]=qbest; //Q
469 for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2;
470 for (Int_t ilab=0;ilab<3;ilab++){
471 milab[ilab] = pos[ip].GetLabel(ilab);
472 milab[ilab+3] = neg[in].GetLabel(ilab);
473 }
474 //
475 CheckLabels2(milab);
476 ratio = (pos[ip].GetQ()-neg[in].GetQ())/(pos[ip].GetQ()+neg[in].GetQ());
477 milab[3]=(((ip<<10) + in)<<10) + idet; // pos|neg|det
478 Int_t info[3] = {pos[ip].GetNd(),neg[in].GetNd(),fNlayer[fModule]};
479
00a7cc50 480 AliITSRecPoint * cl2;
481 if(clusters){
482 cl2 = new (cl[ncl]) AliITSRecPoint(fModule,fDetTypeRec->GetITSgeom(),milab,lp,info);
483 cl2->SetChargeRatio(ratio);
484 cl2->SetType(5);
485 pairs[ip][in] = 5;
486 if ((pos[ip].GetNd()+neg[in].GetNd())>6){ //multi cluster
487 cl2->SetType(6);
488 pairs[ip][in] = 6;
489 }
490 }
04366a57 491 else{
00a7cc50 492 cl2 = new AliITSRecPoint(fModule,fDetTypeRec->GetITSgeom(),milab,lp,info);
493 cl2->SetChargeRatio(ratio);
494 cl2->SetType(5);
495 pairs[ip][in] = 5;
496 if ((pos[ip].GetNd()+neg[in].GetNd())>6){ //multi cluster
497 cl2->SetType(6);
498 pairs[ip][in] = 6;
499 }
500
501 fDetTypeRec->AddRecPoint(*cl2);
04366a57 502 }
503 ncl++;
04366a57 504 }
00a7cc50 505
04366a57 506 //
507 // add second pair
508
00a7cc50 509 // if (!(cused1[ip2] || cused2[in])){ //
04366a57 510 if (pairs[ip2][in]==100){
511 Float_t yp=pos[ip2].GetY()*fYpitchSSD;
512 Float_t yn=neg[in].GetY()*fYpitchSSD;
513 Float_t zt=(2*fHlSSD*tanp + yp - yn)/(tann+tanp);
514 Float_t yt=yn + tann*zt;
515 zt-=fHlSSD; yt-=fHwSSD;
516 ybest =yt; zbest=zt;
517 qbest =pos[ip2].GetQ();
518 lp[0]=-(-ybest+fYshift[fModule]);
519 lp[1]= -zbest+fZshift[fModule];
520 lp[2]=0.0025*0.0025; //SigmaY2
521 lp[3]=0.110*0.110; //SigmaZ2
522
523 lp[4]=qbest; //Q
524 for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2;
525 for (Int_t ilab=0;ilab<3;ilab++){
526 milab[ilab] = pos[ip2].GetLabel(ilab);
527 milab[ilab+3] = neg[in].GetLabel(ilab);
528 }
529 //
530 CheckLabels2(milab);
531 ratio = (pos[ip2].GetQ()-neg[in].GetQ())/(pos[ip2].GetQ()+neg[in].GetQ());
532 milab[3]=(((ip2<<10) + in)<<10) + idet; // pos|neg|det
533 Int_t info[3] = {pos[ip2].GetNd(),neg[in].GetNd(),fNlayer[fModule]};
00a7cc50 534
535 AliITSRecPoint * cl2;
536 if(clusters){
537 cl2 = new (cl[ncl]) AliITSRecPoint(fModule,fDetTypeRec->GetITSgeom(),milab,lp,info);
538 cl2->SetChargeRatio(ratio);
539 cl2->SetType(5);
540 pairs[ip2][in] =5;
541 if ((pos[ip2].GetNd()+neg[in].GetNd())>6){ //multi cluster
542 cl2->SetType(6);
543 pairs[ip2][in] =6;
544 }
545 }
04366a57 546 else{
00a7cc50 547 cl2 = new AliITSRecPoint(fModule,fDetTypeRec->GetITSgeom(),milab,lp,info);
548 cl2->SetChargeRatio(ratio);
549 cl2->SetType(5);
550 pairs[ip2][in] =5;
551 if ((pos[ip2].GetNd()+neg[in].GetNd())>6){ //multi cluster
552 cl2->SetType(6);
553 pairs[ip2][in] =6;
554 }
555
556 fDetTypeRec->AddRecPoint(*cl2);
04366a57 557 }
558 ncl++;
04366a57 559 }
560 cused1[ip]++;
561 cused1[ip2]++;
562 cused2[in]++;
563 }
564 }
565 }
00a7cc50 566
04366a57 567 //
568 for (Int_t jn=0;jn<nn;jn++){
569 if (cused2[jn]) continue;
570 Float_t ybest=1000,zbest=1000,qbest=0;
571 // select "silber" cluster
572 if ( cpositive[jn]==1 && cnegative[positivepair[10*jn]]==2){
573 Int_t ip = positivepair[10*jn];
574 Int_t jn2 = negativepair[10*ip];
575 if (jn2==jn) jn2 = negativepair[10*ip+1];
576 Float_t pcharge = neg[jn].GetQ()+neg[jn2].GetQ();
577 //
578 if (TMath::Abs(pcharge-pos[ip].GetQ())<10){
579 //
580 // add first pair
581 // if (!(cused1[ip]||cused2[jn])){
582 if (pairs[ip][jn]==100){
583 Float_t yn=neg[jn].GetY()*fYpitchSSD;
584 Float_t yp=pos[ip].GetY()*fYpitchSSD;
585 Float_t zt=(2*fHlSSD*tanp + yp - yn)/(tann+tanp);
586 Float_t yt=yn + tann*zt;
587 zt-=fHlSSD; yt-=fHwSSD;
588 ybest =yt; zbest=zt;
589 qbest =neg[jn].GetQ();
590 lp[0]=-(-ybest+fYshift[fModule]);
591 lp[1]= -zbest+fZshift[fModule];
592 lp[2]=0.0025*0.0025; //SigmaY2
593 lp[3]=0.110*0.110; //SigmaZ2
594
595 lp[4]=qbest; //Q
596 for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2;
597 for (Int_t ilab=0;ilab<3;ilab++){
598 milab[ilab] = pos[ip].GetLabel(ilab);
599 milab[ilab+3] = neg[jn].GetLabel(ilab);
600 }
601 //
602 CheckLabels2(milab);
603 ratio = (pos[ip].GetQ()-neg[jn].GetQ())/(pos[ip].GetQ()+neg[jn].GetQ());
604 milab[3]=(((ip<<10) + jn)<<10) + idet; // pos|neg|det
605 Int_t info[3] = {pos[ip].GetNd(),neg[jn].GetNd(),fNlayer[fModule]};
606
00a7cc50 607 AliITSRecPoint * cl2;
608 if(clusters){
609 cl2 = new (cl[ncl]) AliITSRecPoint(fModule,fDetTypeRec->GetITSgeom(),milab,lp,info);
610 cl2->SetChargeRatio(ratio);
611 cl2->SetType(7);
612 pairs[ip][jn] =7;
613 if ((pos[ip].GetNd()+neg[jn].GetNd())>6){ //multi cluster
614 cl2->SetType(8);
615 pairs[ip][jn]=8;
616 }
617
618 }
04366a57 619 else{
00a7cc50 620 cl2 = new AliITSRecPoint(fModule,fDetTypeRec->GetITSgeom(),milab,lp,info);
621 cl2->SetChargeRatio(ratio);
622 cl2->SetType(7);
623 pairs[ip][jn] =7;
624 if ((pos[ip].GetNd()+neg[jn].GetNd())>6){ //multi cluster
625 cl2->SetType(8);
626 pairs[ip][jn]=8;
627 }
628
629 fDetTypeRec->AddRecPoint(*cl2);
04366a57 630 }
631 ncl++;
04366a57 632 }
633 //
634 // add second pair
635 // if (!(cused1[ip]||cused2[jn2])){
636 if (pairs[ip][jn2]==100){
637 Float_t yn=neg[jn2].GetY()*fYpitchSSD;
638 Double_t yp=pos[ip].GetY()*fYpitchSSD;
639 Double_t zt=(2*fHlSSD*tanp + yp - yn)/(tann+tanp);
640 Double_t yt=yn + tann*zt;
641 zt-=fHlSSD; yt-=fHwSSD;
642 ybest =yt; zbest=zt;
643 qbest =neg[jn2].GetQ();
644 lp[0]=-(-ybest+fYshift[fModule]);
645 lp[1]= -zbest+fZshift[fModule];
646 lp[2]=0.0025*0.0025; //SigmaY2
647 lp[3]=0.110*0.110; //SigmaZ2
648
649 lp[4]=qbest; //Q
650 for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2;
651 for (Int_t ilab=0;ilab<3;ilab++){
652 milab[ilab] = pos[ip].GetLabel(ilab);
653 milab[ilab+3] = neg[jn2].GetLabel(ilab);
654 }
655 //
656 CheckLabels2(milab);
657 ratio = (pos[ip].GetQ()-neg[jn2].GetQ())/(pos[ip].GetQ()+neg[jn2].GetQ());
658 milab[3]=(((ip<<10) + jn2)<<10) + idet; // pos|neg|det
659 Int_t info[3] = {pos[ip].GetNd(),neg[jn2].GetNd(),fNlayer[fModule]};
00a7cc50 660 AliITSRecPoint * cl2;
661 if(clusters){
662 cl2 = new (cl[ncl]) AliITSRecPoint(fModule,fDetTypeRec->GetITSgeom(),milab,lp,info);
663 cl2->SetChargeRatio(ratio);
664 pairs[ip][jn2]=7;
665 cl2->SetType(7);
666 if ((pos[ip].GetNd()+neg[jn2].GetNd())>6){ //multi cluster
667 cl2->SetType(8);
668 pairs[ip][jn2]=8;
669 }
670
671 }
04366a57 672 else{
00a7cc50 673 cl2 = new AliITSRecPoint(fModule,fDetTypeRec->GetITSgeom(),milab,lp,info);
674 cl2->SetChargeRatio(ratio);
675 pairs[ip][jn2]=7;
676 cl2->SetType(7);
677 if ((pos[ip].GetNd()+neg[jn2].GetNd())>6){ //multi cluster
678 cl2->SetType(8);
679 pairs[ip][jn2]=8;
680 }
681
682 fDetTypeRec->AddRecPoint(*cl2);
04366a57 683 }
684
685 ncl++;
04366a57 686 }
687 cused1[ip]++;
688 cused2[jn]++;
689 cused2[jn2]++;
690 }
691 }
692 }
693
694 for (Int_t ip=0;ip<np;ip++){
695 Float_t ybest=1000,zbest=1000,qbest=0;
696 //
697 // 2x2 clusters
698 //
699 if ( (cnegative[ip]<5) && cpositive[negativepair[10*ip]]<5){
700 Float_t minchargediff =4.;
701 Int_t j=-1;
702 for (Int_t di=0;di<cnegative[ip];di++){
703 Int_t jc = negativepair[ip*10+di];
704 Float_t chargedif = pos[ip].GetQ()-neg[jc].GetQ();
705 if (TMath::Abs(chargedif)<minchargediff){
706 j =jc;
707 minchargediff = TMath::Abs(chargedif);
708 }
709 }
710 if (j<0) continue; // not proper cluster
711 Int_t count =0;
712 for (Int_t di=0;di<cnegative[ip];di++){
713 Int_t jc = negativepair[ip*10+di];
714 Float_t chargedif = pos[ip].GetQ()-neg[jc].GetQ();
715 if (TMath::Abs(chargedif)<minchargediff+3.) count++;
716 }
717 if (count>1) continue; // more than one "proper" cluster for positive
718 //
719 count =0;
720 for (Int_t dj=0;dj<cpositive[j];dj++){
721 Int_t ic = positivepair[j*10+dj];
722 Float_t chargedif = pos[ic].GetQ()-neg[j].GetQ();
723 if (TMath::Abs(chargedif)<minchargediff+3.) count++;
724 }
725 if (count>1) continue; // more than one "proper" cluster for negative
726
727 Int_t jp = 0;
728
729 count =0;
730 for (Int_t dj=0;dj<cnegative[jp];dj++){
731 Int_t ic = positivepair[jp*10+dj];
732 Float_t chargedif = pos[ic].GetQ()-neg[jp].GetQ();
733 if (TMath::Abs(chargedif)<minchargediff+4.) count++;
734 }
735 if (count>1) continue;
736 if (pairs[ip][j]<100) continue;
737 //
738 //almost gold clusters
739 Float_t yp=pos[ip].GetY()*fYpitchSSD;
740 Float_t yn=neg[j].GetY()*fYpitchSSD;
741 Float_t zt=(2*fHlSSD*tanp + yp - yn)/(tann+tanp);
742 Float_t yt=yn + tann*zt;
743 zt-=fHlSSD; yt-=fHwSSD;
744 ybest=yt; zbest=zt;
745 qbest=0.5*(pos[ip].GetQ()+neg[j].GetQ());
746 lp[0]=-(-ybest+fYshift[fModule]);
747 lp[1]= -zbest+fZshift[fModule];
748 lp[2]=0.0025*0.0025; //SigmaY2
749 lp[3]=0.110*0.110; //SigmaZ2
750 lp[4]=qbest; //Q
751 for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2;
752 for (Int_t ilab=0;ilab<3;ilab++){
753 milab[ilab] = pos[ip].GetLabel(ilab);
754 milab[ilab+3] = neg[j].GetLabel(ilab);
755 }
756 //
757 CheckLabels2(milab);
758 ratio = (pos[ip].GetQ()-neg[j].GetQ())/(pos[ip].GetQ()+neg[j].GetQ());
759 milab[3]=(((ip<<10) + j)<<10) + idet; // pos|neg|det
760 Int_t info[3] = {pos[ip].GetNd(),neg[j].GetNd(),fNlayer[fModule]};
00a7cc50 761 AliITSRecPoint * cl2;
762 if(clusters){
763 cl2 = new (cl[ncl]) AliITSRecPoint(fModule,fDetTypeRec->GetITSgeom(),milab,lp,info);
764 cl2->SetChargeRatio(ratio);
765 cl2->SetType(10);
766 pairs[ip][j]=10;
767 if ((pos[ip].GetNd()+neg[j].GetNd())>6){ //multi cluster
768 cl2->SetType(11);
769 pairs[ip][j]=11;
770 }
771 cused1[ip]++;
772 cused2[j]++;
04366a57 773 }
00a7cc50 774 else{
775 cl2 = new AliITSRecPoint(fModule,fDetTypeRec->GetITSgeom(),milab,lp,info);
776 cl2->SetChargeRatio(ratio);
777 cl2->SetType(10);
778 pairs[ip][j]=10;
779 if ((pos[ip].GetNd()+neg[j].GetNd())>6){ //multi cluster
780 cl2->SetType(11);
781 pairs[ip][j]=11;
782 }
783 cused1[ip]++;
784 cused2[j]++;
785
786 fDetTypeRec->AddRecPoint(*cl2);
787 }
04366a57 788 ncl++;
04366a57 789 }
790
791 }
792
793 //
794 for (Int_t i=0; i<np; i++) {
795 Float_t ybest=1000,zbest=1000,qbest=0;
796 Float_t yp=pos[i].GetY()*fYpitchSSD;
797 if (pos[i].GetQ()<3) continue;
798 for (Int_t j=0; j<nn; j++) {
799 // for (Int_t di = 0;di<cpositive[i];di++){
800 // Int_t j = negativepair[10*i+di];
801 if (neg[j].GetQ()<3) continue;
802 if (cused2[j]||cused1[i]) continue;
803 if (pairs[i][j]>0 &&pairs[i][j]<100) continue;
804 ratio = (pos[i].GetQ()-neg[j].GetQ())/(pos[i].GetQ()+neg[j].GetQ());
805 Float_t yn=neg[j].GetY()*fYpitchSSD;
806 Float_t zt=(2*fHlSSD*tanp + yp - yn)/(tann+tanp);
807 Float_t yt=yn + tann*zt;
808 zt-=fHlSSD; yt-=fHwSSD;
809 if (TMath::Abs(yt)<fHwSSD+0.01)
810 if (TMath::Abs(zt)<fHlSSD+0.01*(neg[j].GetNd()+pos[i].GetNd())) {
811 ybest=yt; zbest=zt;
812 qbest=0.5*(pos[i].GetQ()+neg[j].GetQ());
813 lp[0]=-(-ybest+fYshift[fModule]);
814 lp[1]= -zbest+fZshift[fModule];
815 lp[2]=0.0025*0.0025; //SigmaY2
816 lp[3]=0.110*0.110; //SigmaZ2
817
818 lp[4]=qbest; //Q
819 for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2;
820 for (Int_t ilab=0;ilab<3;ilab++){
821 milab[ilab] = pos[i].GetLabel(ilab);
822 milab[ilab+3] = neg[j].GetLabel(ilab);
823 }
824 //
825 CheckLabels2(milab);
826 milab[3]=(((i<<10) + j)<<10) + idet; // pos|neg|det
827 Int_t info[3] = {pos[i].GetNd(),neg[j].GetNd(),fNlayer[fModule]};
00a7cc50 828 AliITSRecPoint * cl2;
829 if(clusters){
830 cl2 = new (cl[ncl]) AliITSRecPoint(fModule,fDetTypeRec->GetITSgeom(),milab,lp,info);
831 cl2->SetChargeRatio(ratio);
832 cl2->SetType(100+cpositive[j]+cnegative[i]);
833 }
04366a57 834 else{
00a7cc50 835 cl2 = new AliITSRecPoint(fModule,fDetTypeRec->GetITSgeom(),milab,lp,info);
836 cl2->SetChargeRatio(ratio);
837 cl2->SetType(100+cpositive[j]+cnegative[i]);
838 fDetTypeRec->AddRecPoint(*cl2);
04366a57 839 }
840 ncl++;
04366a57 841 //cl2->SetType(0);
842 /*
843 if (pairs[i][j]<100){
844 printf("problem:- %d\n", pairs[i][j]);
845 }
846 if (cnegative[i]<2&&cpositive[j]<2){
847 printf("problem:- %d\n", pairs[i][j]);
848 }
849 */
850 }
851 }
852 }
853
854// for (Int_t i=0; i<1000; i++) delete [] pairs[i];
855// delete [] pairs;
856
857}
858
859
860