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