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