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