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