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