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