]> git.uio.no Git - u/mrichter/AliRoot.git/blame_incremental - ITS/AliITSClusterFinderV2SSD.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame_incremental
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (non-incremental) | history | HEAD
Incremented class version number for SSD calibration objects (Panos)
[u/mrichter/AliRoot.git] / ITS / AliITSClusterFinderV2SSD.cxx
... / ...
CommitLineData
1/**************************************************************************
2 * Copyright(c) 2007-2009, 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/* $Id$ */
17
18////////////////////////////////////////////////////////////////////////////
19// Implementation of the ITS clusterer V2 class //
20// //
21// Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch //
22// Revised: Enrico Fragiacomo, enrico.fragiacomo@ts.infn.it //
23// Revised 23/06/08: Marco Bregant
24// //
25///////////////////////////////////////////////////////////////////////////
26
27#include <Riostream.h>
28
29
30#include "AliITSClusterFinderV2SSD.h"
31#include "AliITSRecPoint.h"
32#include "AliITSgeomTGeo.h"
33#include "AliITSDetTypeRec.h"
34#include "AliRawReader.h"
35#include "AliITSRawStreamSSD.h"
36#include <TClonesArray.h>
37#include "AliITSdigitSSD.h"
38#include "AliITSReconstructor.h"
39#include "AliITSCalibrationSSD.h"
40
41Short_t *AliITSClusterFinderV2SSD::fgPairs = 0x0;
42Int_t AliITSClusterFinderV2SSD::fgPairsSize = 0;
43
44ClassImp(AliITSClusterFinderV2SSD)
45
46
47AliITSClusterFinderV2SSD::AliITSClusterFinderV2SSD(AliITSDetTypeRec* dettyp):AliITSClusterFinderV2(dettyp),
48fLastSSD1(AliITSgeomTGeo::GetModuleIndex(6,1,1)-1),
49fYpitchSSD(0.0095),
50fHwSSD(3.65),
51fHlSSD(2.00),
52fTanP(0.0275),
53fTanN(0.0075)
54{
55
56 //Default constructor
57
58}
59
60//______________________________________________________________________
61AliITSClusterFinderV2SSD::AliITSClusterFinderV2SSD(const AliITSClusterFinderV2SSD &cf) : AliITSClusterFinderV2(cf), fLastSSD1(cf.fLastSSD1),
62fYpitchSSD(cf.fYpitchSSD),
63fHwSSD(cf.fHwSSD),
64fHlSSD(cf.fHlSSD),
65fTanP(cf.fTanP),
66fTanN(cf.fTanN)
67{
68 // Copy constructor
69}
70
71//______________________________________________________________________
72AliITSClusterFinderV2SSD& AliITSClusterFinderV2SSD::operator=(const AliITSClusterFinderV2SSD& cf ){
73 // Assignment operator
74
75 this->~AliITSClusterFinderV2SSD();
76 new(this) AliITSClusterFinderV2SSD(cf);
77 return *this;
78}
79
80
81void AliITSClusterFinderV2SSD::FindRawClusters(Int_t mod){
82
83 //Find clusters V2
84 SetModule(mod);
85 FindClustersSSD(fDigits);
86
87}
88
89void AliITSClusterFinderV2SSD::FindClustersSSD(TClonesArray *alldigits) {
90 //------------------------------------------------------------
91 // Actual SSD cluster finder
92 //------------------------------------------------------------
93
94 static AliITSRecoParam *repa = NULL;
95
96
97 if(!repa){
98 repa = (AliITSRecoParam*) AliITSReconstructor::GetRecoParam();
99 if(!repa){
100 repa = AliITSRecoParam::GetHighFluxParam();
101 AliWarning("Using default AliITSRecoParam class");
102 }
103 }
104
105 AliITSCalibrationSSD* cal = (AliITSCalibrationSSD*)GetResp(fModule);
106 Float_t gain=0;
107
108 Int_t smaxall=alldigits->GetEntriesFast();
109 if (smaxall==0) return;
110 // TObjArray *digits = new TObjArray;
111 TObjArray digits;
112 for (Int_t i=0;i<smaxall; i++){
113 AliITSdigitSSD *d=(AliITSdigitSSD*)alldigits->UncheckedAt(i);
114
115 if(d->IsSideP()) gain = cal->GetGainP(d->GetStripNumber());
116 else gain = cal->GetGainN(d->GetStripNumber());
117
118 Float_t q=gain*d->GetSignal(); // calibration brings mip peaks around 120 (in ADC units)
119 q=cal->ADCToKeV(q); // converts the charge in KeV from ADC units
120 //Float_t q=d->GetSignal()/4.29;// temp. fix (for PID purposed - normalis. to be checked)
121 d->SetSignal(Int_t(q));
122
123 if (d->GetSignal()<3) continue;
124 digits.AddLast(d);
125 }
126 Int_t smax = digits.GetEntriesFast();
127 if (smax==0) return;
128
129 const Int_t kMax=1000;
130 Int_t np=0, nn=0;
131 Ali1Dcluster pos[kMax], neg[kMax];
132 Float_t y=0., q=0., qmax=0.;
133 Int_t lab[4]={-2,-2,-2,-2};
134
135 AliITSdigitSSD *d=(AliITSdigitSSD*)digits.UncheckedAt(0);
136 q += d->GetSignal();
137 y += d->GetCoord2()*d->GetSignal();
138 qmax=d->GetSignal();
139 lab[0]=d->GetTrack(0); lab[1]=d->GetTrack(1); lab[2]=d->GetTrack(2);
140 Int_t curr=d->GetCoord2();
141 Int_t flag=d->GetCoord1();
142 Int_t *n=&nn;
143 Ali1Dcluster *c=neg;
144 Int_t nd=1;
145 Int_t milab[10];
146 for (Int_t ilab=0;ilab<10;ilab++){
147 milab[ilab]=-2;
148 }
149 milab[0]=d->GetTrack(0); milab[1]=d->GetTrack(1); milab[2]=d->GetTrack(2);
150
151 for (Int_t s=1; s<smax; s++) {
152 d=(AliITSdigitSSD*)digits.UncheckedAt(s);
153 Int_t strip=d->GetCoord2();
154 if ((strip-curr) > 1 || flag!=d->GetCoord1()) {
155 c[*n].SetY(y/q);
156 c[*n].SetQ(q);
157 c[*n].SetNd(nd);
158 CheckLabels2(milab);
159 c[*n].SetLabels(milab);
160
161 if(repa->GetUseUnfoldingInClusterFinderSSD()==kTRUE) {
162
163 //Split suspiciously big cluster
164 if (nd>4&&nd<25) {
165 c[*n].SetY(y/q-0.25*nd);
166 c[*n].SetQ(0.5*q);
167 (*n)++;
168 if (*n==kMax) {
169 Error("FindClustersSSD","Too many 1D clusters !");
170 return;
171 }
172 c[*n].SetY(y/q+0.25*nd);
173 c[*n].SetQ(0.5*q);
174 c[*n].SetNd(nd);
175 c[*n].SetLabels(milab);
176 }
177
178 } // unfolding is on
179
180 (*n)++;
181 if (*n==kMax) {
182 Error("FindClustersSSD","Too many 1D clusters !");
183 return;
184 }
185 y=q=qmax=0.;
186 nd=0;
187 lab[0]=lab[1]=lab[2]=-2;
188 //
189 for (Int_t ilab=0;ilab<10;ilab++){
190 milab[ilab]=-2;
191 }
192 //
193 if (flag!=d->GetCoord1()) { n=&np; c=pos; }
194 }
195 flag=d->GetCoord1();
196 q += d->GetSignal();
197 y += d->GetCoord2()*d->GetSignal();
198 nd++;
199 if (d->GetSignal()>qmax) {
200 qmax=d->GetSignal();
201 lab[0]=d->GetTrack(0); lab[1]=d->GetTrack(1); lab[2]=d->GetTrack(2);
202 }
203 for (Int_t ilab=0;ilab<10;ilab++) {
204 if (d->GetTrack(ilab)>=0) AddLabel(milab, (d->GetTrack(ilab)));
205 }
206 curr=strip;
207 }
208 c[*n].SetY(y/q);
209 c[*n].SetQ(q);
210 c[*n].SetNd(nd);
211 c[*n].SetLabels(lab);
212
213 if(repa->GetUseUnfoldingInClusterFinderSSD()==kTRUE) {
214
215 //Split suspiciously big cluster
216 if (nd>4 && nd<25) {
217 c[*n].SetY(y/q-0.25*nd);
218 c[*n].SetQ(0.5*q);
219 (*n)++;
220 if (*n==kMax) {
221 Error("FindClustersSSD","Too many 1D clusters !");
222 return;
223 }
224 c[*n].SetY(y/q+0.25*nd);
225 c[*n].SetQ(0.5*q);
226 c[*n].SetNd(nd);
227 c[*n].SetLabels(lab);
228 }
229 } // unfolding is on
230
231 (*n)++;
232 if (*n==kMax) {
233 Error("FindClustersSSD","Too many 1D clusters !");
234 return;
235 }
236
237 FindClustersSSD(neg, nn, pos, np);
238}
239
240
241void AliITSClusterFinderV2SSD::RawdataToClusters(AliRawReader* rawReader,TClonesArray** clusters){
242
243 //------------------------------------------------------------
244 // This function creates ITS clusters from raw data
245 //------------------------------------------------------------
246 rawReader->Reset();
247 AliITSRawStreamSSD inputSSD(rawReader);
248 FindClustersSSD(&inputSSD,clusters);
249
250}
251
252void AliITSClusterFinderV2SSD::FindClustersSSD(AliITSRawStreamSSD* input,
253 TClonesArray** clusters)
254{
255 //------------------------------------------------------------
256 // Actual SSD cluster finder for raw data
257 //------------------------------------------------------------
258
259 static AliITSRecoParam *repa = NULL;
260 if(!repa){
261 repa = (AliITSRecoParam*) AliITSReconstructor::GetRecoParam();
262 if(!repa){
263 repa = AliITSRecoParam::GetHighFluxParam();
264 AliWarning("Using default AliITSRecoParam class");
265 }
266 }
267
268 Int_t nClustersSSD = 0;
269 const Int_t kMax = 1000;
270 Ali1Dcluster clusters1D[2][kMax];
271 Int_t nClusters[2] = {0, 0};
272 Int_t lab[3]={-2,-2,-2};
273 Float_t q = 0.;
274 Float_t y = 0.;
275 Int_t nDigits = 0;
276 Float_t gain=0;
277 Float_t noise=0.;
278 // Float_t pedestal=0.;
279 Float_t oldnoise=0.;
280 AliITSCalibrationSSD* cal=NULL;
281
282 Int_t matrix[12][1536];
283 Int_t iddl=-1;
284 Int_t iad=-1;
285 Int_t oddl = -1;
286 Int_t oad = -1;
287 Int_t oadc = -1;
288 Int_t ostrip = -1;
289 Int_t osignal = 65535;
290 Int_t n=0;
291 Bool_t next=0;
292
293 // read raw data input stream
294 while (kTRUE) {
295
296 // reset signal matrix
297 for(Int_t i=0; i<12; i++) { for(Int_t j=0; j<1536; j++) { matrix[i][j] = 65535;} }
298
299 if(osignal!=65535) {
300 n++;
301 matrix[oadc][ostrip] = osignal; // recover data from previous occurence of input->Next()
302 }
303
304 // buffer data for ddl=iddl and ad=iad
305 while(kTRUE) {
306
307 next = input->Next();
308 if((!next)&&(input->flag)) continue;
309 Int_t ddl=input->GetDDL();
310 Int_t ad=input->GetAD();
311 Int_t adc = input->GetADC(); adc = (adc<6)? adc : adc - 2;
312 Int_t strip = input->GetStrip();
313 if(input->GetSideFlag()) strip=1535-strip;
314 Int_t signal = input->GetSignal();
315
316 if((ddl==iddl)&&(ad==iad)) {n++; matrix[adc][strip] = signal;}
317 else {oddl=iddl; oad=iad; oadc = adc; ostrip = strip; osignal=signal; iddl=ddl; iad=ad; break;}
318
319 if(!next) {oddl=iddl; oad=iad; oadc = adc; ostrip = strip; osignal=signal; iddl=ddl; iad=ad; break;}
320 //break;
321 }
322
323 // No SSD data
324 if(!next && oddl<0) break;
325
326 if(n==0) continue; // first occurence
327 n=0; //osignal=0;
328
329 // fill 1Dclusters
330 for(Int_t iadc=0; iadc<12; iadc++) { // loop over ADC index for ddl=oddl and ad=oad
331
332 Int_t iimod = (oad - 1) * 12 + iadc;
333 Int_t iModule = AliITSRawStreamSSD::GetModuleNumber(oddl,iimod);
334 if(iModule==-1) continue;
335 cal = (AliITSCalibrationSSD*)GetResp(iModule);
336
337 Bool_t first = 0;
338
339 /*
340 for(Int_t istrip=0; istrip<768; istrip++) { // P-side
341 Int_t signal = matrix[iadc][istrip];
342 pedestal = cal->GetPedestalP(istrip);
343 matrix[iadc][istrip]=signal-(Int_t)pedestal;
344 }
345 */
346
347 /*
348 Float_t cmode=0;
349 for(Int_t l=0; l<6; l++) {
350 cmode=0;
351 for(Int_t n=20; n<108; n++) cmode+=matrix[iadc][l*128+n];
352 cmode/=88.;
353 for(Int_t n=0; n<128; n++) matrix[iadc][l*128+n]-=(Int_t)cmode;
354
355 }
356 */
357
358 Int_t istrip=0;
359 for(istrip=0; istrip<768; istrip++) { // P-side
360
361 Int_t signal = TMath::Abs(matrix[iadc][istrip]);
362
363 oldnoise = noise;
364 noise = cal->GetNoiseP(istrip); if(noise<1.) signal = 65535;
365 if(signal<3*noise) signal = 65535; // in case ZS was not done in hw do it now
366
367 // if(cal->IsPChannelBad(istrip)) signal=0;
368
369 if (signal!=65535) {
370 gain = cal->GetGainP(istrip);
371 signal = (Int_t) ( signal * gain ); // signal is corrected for gain
372 signal = (Int_t) cal->ADCToKeV( signal ); // signal is converted in KeV
373
374 q += signal; // add digit to current cluster
375 y += istrip * signal;
376 nDigits++;
377 first=1;
378 }
379
380 else if(first) {
381
382 if ( ( (nDigits==1) && ( (q==0) || (q>5*oldnoise)) ) || (nDigits>1) ) {
383
384 Ali1Dcluster& cluster = clusters1D[0][nClusters[0]++];
385
386 if(q!=0) cluster.SetY(y/q);
387 else cluster.SetY(istrip-1);
388
389 cluster.SetQ(q);
390 cluster.SetNd(nDigits);
391 cluster.SetLabels(lab);
392
393 if(repa->GetUseUnfoldingInClusterFinderSSD()==kTRUE) {
394
395 //Split suspiciously big cluster
396 if (nDigits > 4&&nDigits < 25) {
397 if(q!=0) cluster.SetY(y/q - 0.25*nDigits);
398 else cluster.SetY(istrip-1 - 0.25*nDigits);
399 cluster.SetQ(0.5*q);
400 if (nClusters[0] == kMax) {
401 Error("FindClustersSSD", "Too many 1D clusters !");
402 return;
403 }
404 Ali1Dcluster& cluster2 = clusters1D[0][nClusters[0]++];
405 if(q!=0) cluster2.SetY(y/q + 0.25*nDigits);
406 else cluster2.SetY(istrip-1 + 0.25*nDigits);
407 cluster2.SetQ(0.5*q);
408 cluster2.SetNd(nDigits);
409 cluster2.SetLabels(lab);
410 }
411 } // unfolding is on
412 }
413
414 y = q = 0.;
415 nDigits = 0;
416 first=0;
417 }
418
419 } // loop over strip on P-side
420
421 // if last strip does have signal
422 if(first) {
423
424 if ( ( (nDigits==1) && ( (q==0) || (q>5*oldnoise)) ) || (nDigits>1) ) {
425
426 Ali1Dcluster& cluster = clusters1D[0][nClusters[0]++];
427
428 if(q!=0) cluster.SetY(y/q);
429 else cluster.SetY(istrip-1);
430
431 cluster.SetQ(q);
432 cluster.SetNd(nDigits);
433 cluster.SetLabels(lab);
434
435 if(repa->GetUseUnfoldingInClusterFinderSSD()==kTRUE) {
436
437 //Split suspiciously big cluster
438 if (nDigits > 4&&nDigits < 25) {
439 if(q!=0) cluster.SetY(y/q - 0.25*nDigits);
440 else cluster.SetY(istrip-1 - 0.25*nDigits);
441 cluster.SetQ(0.5*q);
442 if (nClusters[0] == kMax) {
443 Error("FindClustersSSD", "Too many 1D clusters !");
444 return;
445 }
446 Ali1Dcluster& cluster2 = clusters1D[0][nClusters[0]++];
447 if(q!=0) cluster2.SetY(y/q + 0.25*nDigits);
448 else cluster2.SetY(istrip-1 + 0.25*nDigits);
449 cluster2.SetQ(0.5*q);
450 cluster2.SetNd(nDigits);
451 cluster2.SetLabels(lab);
452 }
453 } // unfolding is on
454
455 }
456 y = q = 0.;
457 nDigits = 0;
458 first=0;
459 }
460
461 /*
462 for(Int_t istrip=768; istrip<1536; istrip++) { // P-side
463 Int_t signal = matrix[iadc][istrip];
464 pedestal = cal->GetPedestalN(1535-istrip);
465 matrix[iadc][istrip]=signal-(Int_t)pedestal;
466 }
467 */
468
469 /*
470 for(Int_t l=6; l<12; l++) {
471 Float_t cmode=0;
472 for(Int_t n=20; n<108; n++) cmode+=matrix[iadc][l*128+n];
473 cmode/=88.;
474 for(Int_t n=0; n<128; n++) matrix[iadc][l*128+n]-=(Int_t)cmode;
475 }
476 */
477
478 oldnoise = 0.;
479 noise = 0.;
480 Int_t strip=0;
481 for(Int_t iistrip=768; iistrip<1536; iistrip++) { // N-side
482
483 Int_t signal = TMath::Abs(matrix[iadc][iistrip]);
484 strip = 1535-iistrip;
485
486 oldnoise = noise;
487 noise = cal->GetNoiseN(strip); if(noise<1.) signal=65535;
488
489 // if(cal->IsNChannelBad(strip)) signal=0;
490
491 if(signal<3*noise) signal = 65535; // in case ZS was not done in hw do it now
492
493 if (signal!=65535) {
494 gain = cal->GetGainN(strip);
495 signal = (Int_t) ( signal * gain); // signal is corrected for gain
496 signal = (Int_t) cal->ADCToKeV( signal ); // signal is converted in KeV
497
498 // add digit to current cluster
499 q += signal;
500 y += strip * signal;
501 nDigits++;
502 first=1;
503 }
504
505 else if(first) {
506
507 if ( ( (nDigits==1) && ( (q==0) || (q>5*oldnoise)) ) || (nDigits>1) ) {
508
509 Ali1Dcluster& cluster = clusters1D[1][nClusters[1]++];
510
511 if(q!=0) cluster.SetY(y/q);
512 else cluster.SetY(strip+1);
513
514 cluster.SetQ(q);
515 cluster.SetNd(nDigits);
516 cluster.SetLabels(lab);
517
518 if(repa->GetUseUnfoldingInClusterFinderSSD()==kTRUE) {
519
520 //Split suspiciously big cluster
521 if (nDigits > 4&&nDigits < 25) {
522 cluster.SetY(y/q - 0.25*nDigits);
523 cluster.SetQ(0.5*q);
524 if (nClusters[1] == kMax) {
525 Error("FindClustersSSD", "Too many 1D clusters !");
526 return;
527 }
528 Ali1Dcluster& cluster2 = clusters1D[1][nClusters[1]++];
529 cluster2.SetY(y/q + 0.25*nDigits);
530 cluster2.SetQ(0.5*q);
531 cluster2.SetNd(nDigits);
532 cluster2.SetLabels(lab);
533 }
534 } // unfolding is on
535 }
536
537 y = q = 0.;
538 nDigits = 0;
539 first=0;
540 }
541
542 } // loop over strips on N-side
543
544 if(first) {
545
546 if ( ( (nDigits==1) && ( (q==0) || (q>5*oldnoise)) ) || (nDigits>1) ) {
547
548 Ali1Dcluster& cluster = clusters1D[1][nClusters[1]++];
549
550 if(q!=0) cluster.SetY(y/q);
551 else cluster.SetY(strip+1);
552
553 cluster.SetQ(q);
554 cluster.SetNd(nDigits);
555 cluster.SetLabels(lab);
556
557 if(repa->GetUseUnfoldingInClusterFinderSSD()==kTRUE) {
558
559 //Split suspiciously big cluster
560 if (nDigits > 4&&nDigits < 25) {
561 if(q!=0) cluster.SetY(y/q - 0.25*nDigits);
562 else cluster.SetY(strip+1 - 0.25*nDigits);
563 cluster.SetQ(0.5*q);
564 if (nClusters[1] == kMax) {
565 Error("FindClustersSSD", "Too many 1D clusters !");
566 return;
567 }
568 Ali1Dcluster& cluster2 = clusters1D[1][nClusters[1]++];
569 if(q!=0) cluster2.SetY(y/q + 0.25*nDigits);
570 else cluster2.SetY(strip+1 + 0.25*nDigits);
571 cluster2.SetQ(0.5*q);
572 cluster2.SetNd(nDigits);
573 cluster2.SetLabels(lab);
574 }
575 } // unfolding is on
576 }
577
578 y = q = 0.;
579 nDigits = 0;
580 first=0;
581 }
582
583 // create recpoints
584 if((nClusters[0])&&(nClusters[1])) {
585
586 clusters[iModule] = new TClonesArray("AliITSRecPoint");
587 fModule = iModule;
588 FindClustersSSD(&clusters1D[0][0], nClusters[0],
589 &clusters1D[1][0], nClusters[1], clusters[iModule]);
590 Int_t nClustersn = clusters[iModule]->GetEntriesFast();
591 nClustersSSD += nClustersn;
592 }
593
594 nClusters[0] = nClusters[1] = 0;
595 y = q = 0.;
596 nDigits = 0;
597
598 } // loop over adc
599
600 if(!next) break;
601 }
602
603 Info("FindClustersSSD", "found clusters in ITS SSD: %d", nClustersSSD);
604}
605
606void AliITSClusterFinderV2SSD::
607FindClustersSSD(Ali1Dcluster* neg, Int_t nn,
608 Ali1Dcluster* pos, Int_t np,
609 TClonesArray *clusters) {
610 //------------------------------------------------------------
611 // Actual SSD cluster finder
612 //------------------------------------------------------------
613
614 const TGeoHMatrix *mT2L=AliITSgeomTGeo::GetTracking2LocalMatrix(fModule);
615
616 TClonesArray &cl=*clusters;
617 //
618 // Float_t tanp=fTanP, tann=fTanN;
619 const Float_t kStartXzero=3.64325;
620 const Float_t kDeltaXzero5or6=0.02239;
621 const Float_t kDeltaZ5to6=7.6/7.0;
622
623 Int_t is6=-1;
624
625 if (fModule>fLastSSD1) is6=1;
626 Int_t idet=fNdet[fModule];
627 Int_t ncl=0;
628
629 //
630 Int_t negativepair[30000];
631 Int_t cnegative[3000];
632 Int_t cused1[3000];
633 Int_t positivepair[30000];
634 Int_t cpositive[3000];
635 Int_t cused2[3000];
636 for (Int_t i=0;i<3000;i++) {cnegative[i]=0; cused1[i]=0;}
637 for (Int_t i=0;i<3000;i++) {cpositive[i]=0; cused2[i]=0;}
638 for (Int_t i=0;i<30000;i++) {negativepair[i]=0; positivepair[i]=0;}
639
640 if ((np*nn) > fgPairsSize) {
641
642 if (fgPairs) delete [] fgPairs;
643 fgPairsSize = 4*np*nn;
644 fgPairs = new Short_t[fgPairsSize];
645 }
646 memset(fgPairs,0,sizeof(Short_t)*np*nn);
647
648 //
649 // find available pairs
650 //
651 for (Int_t i=0; i<np; i++) {
652 Float_t yp=pos[i].GetY();
653 if ( (pos[i].GetQ()>0) && (pos[i].GetQ()<3) ) continue;
654 for (Int_t j=0; j<nn; j++) {
655 if ( (neg[j].GetQ()>0) && (neg[j].GetQ()<3) ) continue;
656 Float_t yn=neg[j].GetY();
657
658 Float_t zt=1.9*(yp-yn)/7.0-is6*kDeltaZ5to6;
659 Float_t yt=is6*(kStartXzero-(285*yp + 1045*yn)/140000.0) + kDeltaXzero5or6;
660
661 if (TMath::Abs(yt)<fHwSSD+0.01)
662 if (TMath::Abs(zt)<fHlSSD+0.01*(neg[j].GetNd()+pos[i].GetNd())) {
663 negativepair[i*10+cnegative[i]] =j; //index
664 positivepair[j*10+cpositive[j]] =i;
665 cnegative[i]++; //counters
666 cpositive[j]++;
667 fgPairs[i*nn+j]=100;
668 }
669 }
670 }
671
672 //
673 // try to recover points out of but close to the module boundaries
674 //
675 for (Int_t i=0; i<np; i++) {
676 Float_t yp=pos[i].GetY();
677 if ( (pos[i].GetQ()>0) && (pos[i].GetQ()<3) ) continue;
678 for (Int_t j=0; j<nn; j++) {
679 if ( (neg[j].GetQ()>0) && (neg[j].GetQ()<3) ) continue;
680 // if both 1Dclusters have an other cross continue
681 if (cpositive[j]&&cnegative[i]) continue;
682 Float_t yn=neg[j].GetY();
683
684 Float_t zt=1.9*(yp-yn)/7.0-is6*kDeltaZ5to6;
685 Float_t yt=is6*(kStartXzero-(285*yp + 1045*yn)/140000.0) + kDeltaXzero5or6;
686
687 if (TMath::Abs(yt)<fHwSSD+0.1)
688 if (TMath::Abs(zt)<fHlSSD+0.15) {
689 // tag 1Dcluster (eventually will produce low quality recpoint)
690 if (cnegative[i]==0) pos[i].SetNd(100); // not available pair
691 if (cpositive[j]==0) neg[j].SetNd(100); // not available pair
692 negativepair[i*10+cnegative[i]] =j; //index
693 positivepair[j*10+cpositive[j]] =i;
694 cnegative[i]++; //counters
695 cpositive[j]++;
696 fgPairs[i*nn+j]=100;
697 }
698 }
699 }
700
701 //
702 Float_t lp[5];
703 Int_t milab[10];
704 Double_t ratio;
705
706
707 static AliITSRecoParam *repa = NULL;
708 if(!repa){
709 repa = (AliITSRecoParam*) AliITSReconstructor::GetRecoParam();
710 if(!repa){
711 repa = AliITSRecoParam::GetHighFluxParam();
712 AliWarning("Using default AliITSRecoParam class");
713 }
714 }
715
716 if(repa->GetUseChargeMatchingInClusterFinderSSD()==kTRUE) {
717
718
719 //
720 // sign gold tracks
721 //
722 for (Int_t ip=0;ip<np;ip++){
723 Float_t ybest=1000,zbest=1000,qbest=0;
724 //
725 // select gold clusters
726 if ( (cnegative[ip]==1) && cpositive[negativepair[10*ip]]==1){
727 Float_t yp=pos[ip].GetY();
728 Int_t j = negativepair[10*ip];
729
730 if( (pos[ip].GetQ()==0) && (neg[j].GetQ() ==0) ) {
731 // both bad, hence continue;
732 // mark both as used (to avoid recover at the end)
733 cused1[ip]++;
734 cused2[j]++;
735 continue;
736 }
737
738 ratio = (pos[ip].GetQ()-neg[j].GetQ())/(pos[ip].GetQ()+neg[j].GetQ());
739
740 // charge matching (note that if posQ or negQ is 0 -> ratio=1 and the following condition is met
741 if (TMath::Abs(ratio)>0.2) continue; // note: 0.2=3xsigma_ratio calculated in cosmics tests
742
743 //
744 Float_t yn=neg[j].GetY();
745
746 Float_t zt=1.9*(yp-yn)/7.0-is6*kDeltaZ5to6;
747 Float_t yt=is6*(kStartXzero-(285*yp + 1045*yn)/140000.0) + kDeltaXzero5or6;
748
749 ybest=yt; zbest=zt;
750
751
752 qbest=0.5*(pos[ip].GetQ()+neg[j].GetQ());
753 if( (pos[ip].GetQ()==0)||(neg[ip].GetQ()==0)) qbest*=2; // in case of bad strips on one side keep all charge from the other one
754
755 {
756 Double_t loc[3]={ybest,0.,zbest},trk[3]={0.,0.,0.};
757 mT2L->MasterToLocal(loc,trk);
758 lp[0]=trk[1];
759 lp[1]=trk[2];
760 }
761 lp[2]=0.0025*0.0025; //SigmaY2
762 lp[3]=0.110*0.110; //SigmaZ2
763
764 lp[4]=qbest; //Q
765 for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2;
766 for (Int_t ilab=0;ilab<3;ilab++){
767 milab[ilab] = pos[ip].GetLabel(ilab);
768 milab[ilab+3] = neg[j].GetLabel(ilab);
769 }
770 //
771 CheckLabels2(milab);
772 milab[3]=(((ip<<10) + j)<<10) + idet; // pos|neg|det
773 Int_t info[3] = {pos[ip].GetNd(),neg[j].GetNd(),fNlayer[fModule]};
774 AliITSRecPoint * cl2;
775
776 if(clusters){ // Note clusters != 0 when method is called for rawdata
777
778
779 cl2 = new (cl[ncl]) AliITSRecPoint(milab,lp,info);
780
781 cl2->SetChargeRatio(ratio);
782 cl2->SetType(1);
783 fgPairs[ip*nn+j]=1;
784
785 if ((pos[ip].GetNd()+neg[j].GetNd())>6){ //multi cluster
786 cl2->SetType(2);
787 fgPairs[ip*nn+j]=2;
788 }
789
790 if(pos[ip].GetQ()==0) cl2->SetType(3);
791 if(neg[ip].GetQ()==0) cl2->SetType(4);
792
793 cused1[ip]++;
794 cused2[j]++;
795
796 }
797 else{ // Note clusters == 0 when method is called for digits
798
799 cl2 = new AliITSRecPoint(milab,lp,info);
800
801 cl2->SetChargeRatio(ratio);
802 cl2->SetType(1);
803 fgPairs[ip*nn+j]=1;
804
805 if ((pos[ip].GetNd()+neg[j].GetNd())>6){ //multi cluster
806 cl2->SetType(2);
807 fgPairs[ip*nn+j]=2;
808 }
809
810 if(pos[ip].GetQ()==0) cl2->SetType(3);
811 if(neg[ip].GetQ()==0) cl2->SetType(4);
812
813 cused1[ip]++;
814 cused2[j]++;
815
816 fDetTypeRec->AddRecPoint(*cl2);
817 }
818 ncl++;
819 }
820 }
821
822 for (Int_t ip=0;ip<np;ip++){
823 Float_t ybest=1000,zbest=1000,qbest=0;
824 //
825 //
826 // select "silber" cluster
827 if ( cnegative[ip]==1 && cpositive[negativepair[10*ip]]==2){
828 Int_t in = negativepair[10*ip];
829 Int_t ip2 = positivepair[10*in];
830 if (ip2==ip) ip2 = positivepair[10*in+1];
831 Float_t pcharge = pos[ip].GetQ()+pos[ip2].GetQ();
832
833 if ( (TMath::Abs(pcharge-neg[in].GetQ())<30) && (pcharge!=0) ) { //
834
835 //
836 // add first pair
837 if ( (fgPairs[ip*nn+in]==100)&&(pos[ip].GetQ() ) ) { //
838
839 Float_t yp=pos[ip].GetY();
840 Float_t yn=neg[in].GetY();
841
842 Float_t zt=1.9*(yp-yn)/7.0-is6*kDeltaZ5to6;
843 Float_t yt=is6*(kStartXzero-(285*yp + 1045*yn)/140000.0) + kDeltaXzero5or6;
844 ybest =yt; zbest=zt;
845
846
847 qbest =pos[ip].GetQ();
848 Double_t loc[3]={ybest,0.,zbest},trk[3]={0.,0.,0.};
849 mT2L->MasterToLocal(loc,trk);
850 lp[0]=trk[1];
851 lp[1]=trk[2];
852
853 lp[2]=0.0025*0.0025; //SigmaY2
854 lp[3]=0.110*0.110; //SigmaZ2
855
856 lp[4]=qbest; //Q
857 for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2;
858 for (Int_t ilab=0;ilab<3;ilab++){
859 milab[ilab] = pos[ip].GetLabel(ilab);
860 milab[ilab+3] = neg[in].GetLabel(ilab);
861 }
862 //
863 CheckLabels2(milab);
864 ratio = (pos[ip].GetQ()-neg[in].GetQ())/(pos[ip].GetQ()+neg[in].GetQ());
865 milab[3]=(((ip<<10) + in)<<10) + idet; // pos|neg|det
866 Int_t info[3] = {pos[ip].GetNd(),neg[in].GetNd(),fNlayer[fModule]};
867
868 AliITSRecPoint * cl2;
869 if(clusters){
870
871 cl2 = new (cl[ncl]) AliITSRecPoint(milab,lp,info);
872 cl2->SetChargeRatio(ratio);
873 cl2->SetType(5);
874 fgPairs[ip*nn+in] = 5;
875 if ((pos[ip].GetNd()+neg[in].GetNd())>6){ //multi cluster
876 cl2->SetType(6);
877 fgPairs[ip*nn+in] = 6;
878 }
879 }
880 else{
881 cl2 = new AliITSRecPoint(milab,lp,info);
882 cl2->SetChargeRatio(ratio);
883 cl2->SetType(5);
884 fgPairs[ip*nn+in] = 5;
885 if ((pos[ip].GetNd()+neg[in].GetNd())>6){ //multi cluster
886 cl2->SetType(6);
887 fgPairs[ip*nn+in] = 6;
888 }
889
890 fDetTypeRec->AddRecPoint(*cl2);
891 }
892 ncl++;
893 }
894
895
896 //
897 // add second pair
898
899 // if (!(cused1[ip2] || cused2[in])){ //
900 if ( (fgPairs[ip2*nn+in]==100) && (pos[ip2].GetQ()) ) {
901
902 Float_t yp=pos[ip2].GetY();
903 Float_t yn=neg[in].GetY();
904
905 Float_t zt=1.9*(yp-yn)/7.0-is6*kDeltaZ5to6;
906 Float_t yt=is6*(kStartXzero-(285*yp + 1045*yn)/140000.0) + kDeltaXzero5or6;
907 ybest =yt; zbest=zt;
908
909 qbest =pos[ip2].GetQ();
910
911 Double_t loc[3]={ybest,0.,zbest},trk[3]={0.,0.,0.};
912 mT2L->MasterToLocal(loc,trk);
913 lp[0]=trk[1];
914 lp[1]=trk[2];
915
916 lp[2]=0.0025*0.0025; //SigmaY2
917 lp[3]=0.110*0.110; //SigmaZ2
918
919 lp[4]=qbest; //Q
920 for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2;
921 for (Int_t ilab=0;ilab<3;ilab++){
922 milab[ilab] = pos[ip2].GetLabel(ilab);
923 milab[ilab+3] = neg[in].GetLabel(ilab);
924 }
925 //
926 CheckLabels2(milab);
927 ratio = (pos[ip2].GetQ()-neg[in].GetQ())/(pos[ip2].GetQ()+neg[in].GetQ());
928 milab[3]=(((ip2<<10) + in)<<10) + idet; // pos|neg|det
929 Int_t info[3] = {pos[ip2].GetNd(),neg[in].GetNd(),fNlayer[fModule]};
930
931 AliITSRecPoint * cl2;
932 if(clusters){
933 cl2 = new (cl[ncl]) AliITSRecPoint(milab,lp,info);
934
935 cl2->SetChargeRatio(ratio);
936 cl2->SetType(5);
937 fgPairs[ip2*nn+in] =5;
938 if ((pos[ip2].GetNd()+neg[in].GetNd())>6){ //multi cluster
939 cl2->SetType(6);
940 fgPairs[ip2*nn+in] =6;
941 }
942 }
943 else{
944 cl2 = new AliITSRecPoint(milab,lp,info);
945 cl2->SetChargeRatio(ratio);
946 cl2->SetType(5);
947 fgPairs[ip2*nn+in] =5;
948 if ((pos[ip2].GetNd()+neg[in].GetNd())>6){ //multi cluster
949 cl2->SetType(6);
950 fgPairs[ip2*nn+in] =6;
951 }
952
953 fDetTypeRec->AddRecPoint(*cl2);
954 }
955 ncl++;
956 }
957
958 cused1[ip]++;
959 cused1[ip2]++;
960 cused2[in]++;
961
962 } // charge matching condition
963
964 } // 2 Pside cross 1 Nside
965 } // loop over Pside clusters
966
967
968
969 //
970 for (Int_t jn=0;jn<nn;jn++){
971 if (cused2[jn]) continue;
972 Float_t ybest=1000,zbest=1000,qbest=0;
973 // select "silber" cluster
974 if ( cpositive[jn]==1 && cnegative[positivepair[10*jn]]==2){
975 Int_t ip = positivepair[10*jn];
976 Int_t jn2 = negativepair[10*ip];
977 if (jn2==jn) jn2 = negativepair[10*ip+1];
978 Float_t pcharge = neg[jn].GetQ()+neg[jn2].GetQ();
979 //
980
981 if ( (TMath::Abs(pcharge-pos[ip].GetQ())<30) && // charge matching
982 (pcharge!=0) ) { // reject combinations of bad strips
983
984 //
985 // add first pair
986 // if (!(cused1[ip]||cused2[jn])){
987 if ( (fgPairs[ip*nn+jn]==100) && (neg[jn].GetQ()) ) { //
988
989 Float_t yn=neg[jn].GetY();
990 Float_t yp=pos[ip].GetY();
991
992 Float_t zt=1.9*(yp-yn)/7.0-is6*kDeltaZ5to6;
993 Float_t yt=is6*(kStartXzero-(285*yp + 1045*yn)/140000.0) + kDeltaXzero5or6;
994
995 ybest =yt; zbest=zt;
996
997
998 qbest =neg[jn].GetQ();
999
1000 {
1001 Double_t loc[3]={ybest,0.,zbest},trk[3]={0.,0.,0.};
1002 mT2L->MasterToLocal(loc,trk);
1003 lp[0]=trk[1];
1004 lp[1]=trk[2];
1005 }
1006 lp[2]=0.0025*0.0025; //SigmaY2
1007 lp[3]=0.110*0.110; //SigmaZ2
1008
1009 lp[4]=qbest; //Q
1010 for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2;
1011 for (Int_t ilab=0;ilab<3;ilab++){
1012 milab[ilab] = pos[ip].GetLabel(ilab);
1013 milab[ilab+3] = neg[jn].GetLabel(ilab);
1014 }
1015 //
1016 CheckLabels2(milab);
1017 ratio = (pos[ip].GetQ()-neg[jn].GetQ())/(pos[ip].GetQ()+neg[jn].GetQ());
1018 milab[3]=(((ip<<10) + jn)<<10) + idet; // pos|neg|det
1019 Int_t info[3] = {pos[ip].GetNd(),neg[jn].GetNd(),fNlayer[fModule]};
1020
1021 AliITSRecPoint * cl2;
1022 if(clusters){
1023 cl2 = new (cl[ncl]) AliITSRecPoint(milab,lp,info);
1024
1025 cl2->SetChargeRatio(ratio);
1026 cl2->SetType(7);
1027 fgPairs[ip*nn+jn] =7;
1028 if ((pos[ip].GetNd()+neg[jn].GetNd())>6){ //multi cluster
1029 cl2->SetType(8);
1030 fgPairs[ip*nn+jn]=8;
1031 }
1032
1033 }
1034 else{
1035 cl2 = new AliITSRecPoint(milab,lp,info);
1036 cl2->SetChargeRatio(ratio);
1037 cl2->SetType(7);
1038 fgPairs[ip*nn+jn] =7;
1039 if ((pos[ip].GetNd()+neg[jn].GetNd())>6){ //multi cluster
1040 cl2->SetType(8);
1041 fgPairs[ip*nn+jn]=8;
1042 }
1043
1044 fDetTypeRec->AddRecPoint(*cl2);
1045 }
1046 ncl++;
1047 }
1048 //
1049 // add second pair
1050 // if (!(cused1[ip]||cused2[jn2])){
1051 if ( (fgPairs[ip*nn+jn2]==100)&&(neg[jn2].GetQ() ) ) { //
1052
1053 Float_t yn=neg[jn2].GetY();
1054 Double_t yp=pos[ip].GetY();
1055 Double_t zt=1.9*(yp-yn)/7.0-is6*kDeltaZ5to6;
1056 Double_t yt=is6*(kStartXzero-(285*yp + 1045*yn)/140000.0) + kDeltaXzero5or6;
1057
1058 ybest =yt; zbest=zt;
1059
1060
1061 qbest =neg[jn2].GetQ();
1062
1063 {
1064 Double_t loc[3]={ybest,0.,zbest},trk[3]={0.,0.,0.};
1065 mT2L->MasterToLocal(loc,trk);
1066 lp[0]=trk[1];
1067 lp[1]=trk[2];
1068 }
1069 lp[2]=0.0025*0.0025; //SigmaY2
1070 lp[3]=0.110*0.110; //SigmaZ2
1071
1072 lp[4]=qbest; //Q
1073 for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2;
1074 for (Int_t ilab=0;ilab<3;ilab++){
1075 milab[ilab] = pos[ip].GetLabel(ilab);
1076 milab[ilab+3] = neg[jn2].GetLabel(ilab);
1077 }
1078 //
1079 CheckLabels2(milab);
1080 ratio = (pos[ip].GetQ()-neg[jn2].GetQ())/(pos[ip].GetQ()+neg[jn2].GetQ());
1081 milab[3]=(((ip<<10) + jn2)<<10) + idet; // pos|neg|det
1082 Int_t info[3] = {pos[ip].GetNd(),neg[jn2].GetNd(),fNlayer[fModule]};
1083 AliITSRecPoint * cl2;
1084 if(clusters){
1085 cl2 = new (cl[ncl]) AliITSRecPoint(milab,lp,info);
1086
1087
1088 cl2->SetChargeRatio(ratio);
1089 fgPairs[ip*nn+jn2]=7;
1090 cl2->SetType(7);
1091 if ((pos[ip].GetNd()+neg[jn2].GetNd())>6){ //multi cluster
1092 cl2->SetType(8);
1093 fgPairs[ip*nn+jn2]=8;
1094 }
1095
1096 }
1097 else{
1098 cl2 = new AliITSRecPoint(milab,lp,info);
1099 cl2->SetChargeRatio(ratio);
1100 fgPairs[ip*nn+jn2]=7;
1101 cl2->SetType(7);
1102 if ((pos[ip].GetNd()+neg[jn2].GetNd())>6){ //multi cluster
1103 cl2->SetType(8);
1104 fgPairs[ip*nn+jn2]=8;
1105 }
1106
1107 fDetTypeRec->AddRecPoint(*cl2);
1108 }
1109
1110 ncl++;
1111 }
1112 cused1[ip]++;
1113 cused2[jn]++;
1114 cused2[jn2]++;
1115
1116 } // charge matching condition
1117
1118 } // 2 Nside cross 1 Pside
1119 } // loop over Pside clusters
1120
1121
1122
1123 for (Int_t ip=0;ip<np;ip++){
1124 Float_t ybest=1000,zbest=1000,qbest=0;
1125 //
1126 // 2x2 clusters
1127 //
1128 if ( (cnegative[ip]<5) && cpositive[negativepair[10*ip]]<5){
1129 Float_t minchargediff =4.;
1130 Int_t j=-1;
1131 for (Int_t di=0;di<cnegative[ip];di++){
1132 Int_t jc = negativepair[ip*10+di];
1133 Float_t chargedif = pos[ip].GetQ()-neg[jc].GetQ();
1134 if (TMath::Abs(chargedif)<minchargediff){
1135 j =jc;
1136 minchargediff = TMath::Abs(chargedif);
1137 }
1138 }
1139 if (j<0) continue; // not proper cluster
1140
1141 Int_t count =0;
1142 for (Int_t di=0;di<cnegative[ip];di++){
1143 Int_t jc = negativepair[ip*10+di];
1144 Float_t chargedif = pos[ip].GetQ()-neg[jc].GetQ();
1145 if (TMath::Abs(chargedif)<minchargediff+3.) count++;
1146 }
1147 if (count>1) continue; // more than one "proper" cluster for positive
1148 //
1149
1150 count =0;
1151 for (Int_t dj=0;dj<cpositive[j];dj++){
1152 Int_t ic = positivepair[j*10+dj];
1153 Float_t chargedif = pos[ic].GetQ()-neg[j].GetQ();
1154 if (TMath::Abs(chargedif)<minchargediff+3.) count++;
1155 }
1156 if (count>1) continue; // more than one "proper" cluster for negative
1157
1158 Int_t jp = 0;
1159
1160 count =0;
1161 for (Int_t dj=0;dj<cnegative[jp];dj++){
1162 Int_t ic = positivepair[jp*10+dj];
1163 Float_t chargedif = pos[ic].GetQ()-neg[jp].GetQ();
1164 if (TMath::Abs(chargedif)<minchargediff+4.) count++;
1165 }
1166 if (count>1) continue;
1167 if (fgPairs[ip*nn+j]<100) continue;
1168 //
1169
1170 //almost gold clusters
1171 Float_t yp=pos[ip].GetY();
1172 Float_t yn=neg[j].GetY();
1173
1174 Float_t zt=1.9*(yp-yn)/7.0-is6*kDeltaZ5to6;
1175 Float_t yt=is6*(kStartXzero-(285*yp + 1045*yn)/140000.0) + kDeltaXzero5or6;
1176
1177 ybest=yt; zbest=zt;
1178 qbest=0.5*(pos[ip].GetQ()+neg[j].GetQ());
1179
1180 {
1181 Double_t loc[3]={ybest,0.,zbest},trk[3]={0.,0.,0.};
1182 mT2L->MasterToLocal(loc,trk);
1183 lp[0]=trk[1];
1184 lp[1]=trk[2];
1185 }
1186 lp[2]=0.0025*0.0025; //SigmaY2
1187 lp[3]=0.110*0.110; //SigmaZ2
1188 lp[4]=qbest; //Q
1189 for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2;
1190 for (Int_t ilab=0;ilab<3;ilab++){
1191 milab[ilab] = pos[ip].GetLabel(ilab);
1192 milab[ilab+3] = neg[j].GetLabel(ilab);
1193 }
1194 //
1195 CheckLabels2(milab);
1196 if ((neg[j].GetQ()==0)&&(pos[ip].GetQ()==0)) continue; // reject crosses of bad strips!!
1197 ratio = (pos[ip].GetQ()-neg[j].GetQ())/(pos[ip].GetQ()+neg[j].GetQ());
1198 milab[3]=(((ip<<10) + j)<<10) + idet; // pos|neg|det
1199 Int_t info[3] = {pos[ip].GetNd(),neg[j].GetNd(),fNlayer[fModule]};
1200 AliITSRecPoint * cl2;
1201 if(clusters){
1202 cl2 = new (cl[ncl]) AliITSRecPoint(milab,lp,info);
1203
1204 cl2->SetChargeRatio(ratio);
1205 cl2->SetType(10);
1206 fgPairs[ip*nn+j]=10;
1207 if ((pos[ip].GetNd()+neg[j].GetNd())>6){ //multi cluster
1208 cl2->SetType(11);
1209 fgPairs[ip*nn+j]=11;
1210 }
1211 cused1[ip]++;
1212 cused2[j]++;
1213 }
1214 else{
1215 cl2 = new AliITSRecPoint(milab,lp,info);
1216 cl2->SetChargeRatio(ratio);
1217 cl2->SetType(10);
1218 fgPairs[ip*nn+j]=10;
1219 if ((pos[ip].GetNd()+neg[j].GetNd())>6){ //multi cluster
1220 cl2->SetType(11);
1221 fgPairs[ip*nn+j]=11;
1222 }
1223 cused1[ip]++;
1224 cused2[j]++;
1225
1226 fDetTypeRec->AddRecPoint(*cl2);
1227 }
1228 ncl++;
1229
1230 } // manyXmany
1231 } // loop over Pside 1Dclusters
1232
1233
1234 } // use charge matching
1235
1236
1237 // recover all the other crosses
1238 //
1239 for (Int_t i=0; i<np; i++) {
1240 Float_t ybest=1000,zbest=1000,qbest=0;
1241 Float_t yp=pos[i].GetY();
1242 if ((pos[i].GetQ()>0)&&(pos[i].GetQ()<3)) continue;
1243 for (Int_t j=0; j<nn; j++) {
1244 // for (Int_t di = 0;di<cpositive[i];di++){
1245 // Int_t j = negativepair[10*i+di];
1246 if ((neg[j].GetQ()>0)&&(neg[j].GetQ()<3)) continue;
1247
1248 if ((neg[j].GetQ()==0)&&(pos[i].GetQ()==0)) continue; // reject crosses of bad strips!!
1249
1250 if (cused2[j]||cused1[i]) continue;
1251 if (fgPairs[i*nn+j]>0 &&fgPairs[i*nn+j]<100) continue;
1252 ratio = (pos[i].GetQ()-neg[j].GetQ())/(pos[i].GetQ()+neg[j].GetQ());
1253 Float_t yn=neg[j].GetY();
1254
1255 Float_t zt=1.9*(yp-yn)/7.0-is6*kDeltaZ5to6;
1256 Float_t yt=is6*(kStartXzero-(285*yp + 1045*yn)/140000.0) + kDeltaXzero5or6;
1257
1258 if (TMath::Abs(yt)<fHwSSD+0.01)
1259 if (TMath::Abs(zt)<fHlSSD+0.01*(neg[j].GetNd()+pos[i].GetNd())) {
1260 ybest=yt; zbest=zt;
1261 qbest=0.5*(pos[i].GetQ()+neg[j].GetQ());
1262
1263 {
1264 Double_t loc[3]={ybest,0.,zbest},trk[3]={0.,0.,0.};
1265 mT2L->MasterToLocal(loc,trk);
1266 lp[0]=trk[1];
1267 lp[1]=trk[2];
1268 }
1269 lp[2]=0.0025*0.0025; //SigmaY2
1270 lp[3]=0.110*0.110; //SigmaZ2
1271
1272 lp[4]=qbest; //Q
1273 for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2;
1274 for (Int_t ilab=0;ilab<3;ilab++){
1275 milab[ilab] = pos[i].GetLabel(ilab);
1276 milab[ilab+3] = neg[j].GetLabel(ilab);
1277 }
1278 //
1279 CheckLabels2(milab);
1280 milab[3]=(((i<<10) + j)<<10) + idet; // pos|neg|det
1281 Int_t info[3] = {pos[i].GetNd(),neg[j].GetNd(),fNlayer[fModule]};
1282 AliITSRecPoint * cl2;
1283 if(clusters){
1284 cl2 = new (cl[ncl]) AliITSRecPoint(milab,lp,info);
1285
1286 cl2->SetChargeRatio(ratio);
1287 cl2->SetType(100+cpositive[j]+cnegative[i]);
1288
1289 if(pos[i].GetQ()==0) cl2->SetType(200+cpositive[j]+cnegative[i]);
1290 if(neg[j].GetQ()==0) cl2->SetType(300+cpositive[j]+cnegative[i]);
1291
1292 }
1293 else{
1294 cl2 = new AliITSRecPoint(milab,lp,info);
1295 cl2->SetChargeRatio(ratio);
1296 cl2->SetType(100+cpositive[j]+cnegative[i]);
1297
1298 if(pos[i].GetQ()==0) cl2->SetType(200+cpositive[j]+cnegative[i]);
1299 if(neg[j].GetQ()==0) cl2->SetType(300+cpositive[j]+cnegative[i]);
1300
1301 fDetTypeRec->AddRecPoint(*cl2);
1302 }
1303 ncl++;
1304 }
1305 }
1306 }
1307
1308}