]> 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
Changing a default argument in the standard constructor (M. Van Leeuwen)
[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 //
d695268b 22// Last revision: 13-05-09 Enrico Fragiacomo //
23// enrico.fragiacomo@ts.infn.it //
04366a57 24// //
25///////////////////////////////////////////////////////////////////////////
26
308b5ea4 27#include <Riostream.h>
d695268b 28#include "AliLog.h"
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;
d695268b 44const Float_t AliITSClusterFinderV2SSD::fgkThreshold = 5.;
45
42ed6062 46const Float_t AliITSClusterFinderV2SSD::fgkCosmic2008StripShifts[16][9] =
47 {{-0.35,-0.35,-0.35,-0.35,-0.35,-0.35,-0.35,-0.35,-0.35}, // DDL 512
48 {-0.35,-0.35,-0.35,-0.35,-0.35,-0.35,-0.35,-0.35,-0.35}, // DDL 513
49 {-0.15,-0.15,-0.15,-0.15,-0.15,-0.15,-0.15,-0.15,-0.15}, // DDL 514
50 {-0.15,-0.15,-0.15,-0.15,-0.15,-0.15,-0.15,-0.15,-0.15}, // DDL 515
51 { 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00}, // DDL 516
52 { 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00}, // DDL 517
53 {-0.15,-0.15,-0.15,-0.15,-0.15,-0.15,-0.15,-0.15,-0.15}, // DDL 518
54 {-0.15,-0.15,-0.15,-0.15,-0.15,-0.15,-0.15,-0.15,-0.15}, // DDL 519
55 {-0.15,-0.15,-0.15,-0.15,-0.15,-0.15,-0.15,-0.25,-0.15}, // DDL 520
56 {-0.15,-0.15,-0.15,-0.15,-0.15,-0.15,-0.15,-0.15,-0.15}, // DDL 521
57 {-0.10,-0.10,-0.10,-0.40,-0.40,-0.40,-0.10,-0.10,-0.45}, // DDL 522
58 {-0.10,-0.10,-0.10,-0.35,-0.35,-0.35,-0.10,-0.35,-0.50}, // DDL 523
59 { 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00}, // DDL 524
60 { 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00}, // DDL 525
61 { 0.35, 0.35, 0.35, 0.35, 0.35, 0.35, 0.35, 0.35, 0.35}, // DDL 526
62 { 0.45, 0.45, 0.45, 0.45, 0.45, 0.45, 0.45, 0.45, 0.45}}; // DDL 527
c157c94e 63
04366a57 64ClassImp(AliITSClusterFinderV2SSD)
65
66
5d2c2f86 67AliITSClusterFinderV2SSD::AliITSClusterFinderV2SSD(AliITSDetTypeRec* dettyp):AliITSClusterFinder(dettyp),
42ed6062 68 fLastSSD1(AliITSgeomTGeo::GetModuleIndex(6,1,1)-1)
d036ccd3 69{
8be4e1b1 70//Default constructor
04366a57 71
04366a57 72}
73
308b5ea4 74//______________________________________________________________________
5d2c2f86 75AliITSClusterFinderV2SSD::AliITSClusterFinderV2SSD(const AliITSClusterFinderV2SSD &cf) : AliITSClusterFinder(cf), fLastSSD1(cf.fLastSSD1)
308b5ea4 76{
77 // Copy constructor
308b5ea4 78}
79
80//______________________________________________________________________
81AliITSClusterFinderV2SSD& AliITSClusterFinderV2SSD::operator=(const AliITSClusterFinderV2SSD& cf ){
82 // Assignment operator
83
84 this->~AliITSClusterFinderV2SSD();
85 new(this) AliITSClusterFinderV2SSD(cf);
86 return *this;
87}
88
04366a57 89
90void AliITSClusterFinderV2SSD::FindRawClusters(Int_t mod){
91
92 //Find clusters V2
93 SetModule(mod);
94 FindClustersSSD(fDigits);
95
96}
97
98void AliITSClusterFinderV2SSD::FindClustersSSD(TClonesArray *alldigits) {
99 //------------------------------------------------------------
100 // Actual SSD cluster finder
101 //------------------------------------------------------------
d695268b 102 Int_t smaxall=alldigits->GetEntriesFast();
103 if (smaxall==0) return;
a86176e3 104
d695268b 105
106 //---------------------------------------
107 // load recoparam and calibration
108 //
109 static AliITSRecoParam *repa = NULL;
a86176e3 110 if(!repa){
111 repa = (AliITSRecoParam*) AliITSReconstructor::GetRecoParam();
112 if(!repa){
ed446fa3 113 repa = AliITSRecoParam::GetHighFluxParam();
a86176e3 114 AliWarning("Using default AliITSRecoParam class");
115 }
116 }
117
3a4139a2 118 AliITSCalibrationSSD* cal = (AliITSCalibrationSSD*)GetResp(fModule);
119 Float_t gain=0;
d695268b 120 Float_t noise=0;
121 //---------------------------------------
3a4139a2 122
d695268b 123
124 //------------------------------------
125 // fill the digits array with zero-suppression condition
126 // Signal is converted in KeV
127 //
f7b30404 128 TObjArray digits;
04366a57 129 for (Int_t i=0;i<smaxall; i++){
130 AliITSdigitSSD *d=(AliITSdigitSSD*)alldigits->UncheckedAt(i);
3a4139a2 131
d695268b 132 if(d->IsSideP()) noise = cal->GetNoiseP(d->GetStripNumber());
133 else noise = cal->GetNoiseN(d->GetStripNumber());
134 if (d->GetSignal()<3.*noise) continue;
135
3a4139a2 136 if(d->IsSideP()) gain = cal->GetGainP(d->GetStripNumber());
137 else gain = cal->GetGainN(d->GetStripNumber());
138
d695268b 139 Float_t q=gain*d->GetSignal(); //
3a4139a2 140 q=cal->ADCToKeV(q); // converts the charge in KeV from ADC units
6490e1c5 141 d->SetSignal(Int_t(q));
3a4139a2 142
f7b30404 143 digits.AddLast(d);
04366a57 144 }
f7b30404 145 Int_t smax = digits.GetEntriesFast();
04366a57 146 if (smax==0) return;
d695268b 147 //------------------------------------
148
04366a57 149
150 const Int_t kMax=1000;
151 Int_t np=0, nn=0;
152 Ali1Dcluster pos[kMax], neg[kMax];
153 Float_t y=0., q=0., qmax=0.;
d695268b 154 Int_t lab[4]={-2,-2,-2,-2};
155 Bool_t flag5 = 0;
04366a57 156
d695268b 157 /*
158 cout<<"-----------------------------"<<endl;
159 cout<<"this is module "<<fModule;
160 cout<<endl;
161 cout<<endl;
162 */
163
164 //--------------------------------------------------------
165 // start 1D-clustering from the first digit in the digits array
166 //
f7b30404 167 AliITSdigitSSD *d=(AliITSdigitSSD*)digits.UncheckedAt(0);
04366a57 168 q += d->GetSignal();
169 y += d->GetCoord2()*d->GetSignal();
170 qmax=d->GetSignal();
171 lab[0]=d->GetTrack(0); lab[1]=d->GetTrack(1); lab[2]=d->GetTrack(2);
d695268b 172
173 if(d->IsSideP()) {
174 noise = cal->GetNoiseP(d->GetStripNumber());
175 gain = cal->GetGainP(d->GetStripNumber());
176 }
177 else {
178 noise = cal->GetNoiseN(d->GetStripNumber());
179 gain = cal->GetGainN(d->GetStripNumber());
180 }
181 noise*=gain;
182 noise=cal->ADCToKeV(noise); // converts noise in KeV from ADC units
183
184 if(qmax>fgkThreshold*noise) flag5=1; // seed for the cluster
185
186 /*
187 cout<<d->GetSignal()<<" "<<noise<<" "<<flag5<<" "<<
188 d->GetCoord1()<<" "<<d->GetCoord2()<<endl;
189 */
190
04366a57 191 Int_t curr=d->GetCoord2();
192 Int_t flag=d->GetCoord1();
d695268b 193
194 // Note: the first side which will be processed is supposed to be the
195 // P-side which is neg
04366a57 196 Int_t *n=&nn;
197 Ali1Dcluster *c=neg;
d695268b 198 if(flag) {n=&np; c=pos;} // in case we have only Nstrips (P was bad!)
199
04366a57 200 Int_t nd=1;
201 Int_t milab[10];
202 for (Int_t ilab=0;ilab<10;ilab++){
203 milab[ilab]=-2;
204 }
205 milab[0]=d->GetTrack(0); milab[1]=d->GetTrack(1); milab[2]=d->GetTrack(2);
206
d695268b 207
208 //----------------------------------------------------------
209 // search for neighboring digits
210 //
04366a57 211 for (Int_t s=1; s<smax; s++) {
f7b30404 212 d=(AliITSdigitSSD*)digits.UncheckedAt(s);
04366a57 213 Int_t strip=d->GetCoord2();
d695268b 214
215 // if digits is not a neighbour or side did not change
216 // and at least one of the previous digits met the seed condition
217 // then creates a new 1D cluster
218 if ( ( ((strip-curr) > 1) || (flag!=d->GetCoord1()) ) ) {
219
220 if(flag5) {
221 //cout<<"here1"<<endl;
04366a57 222 c[*n].SetY(y/q);
223 c[*n].SetQ(q);
224 c[*n].SetNd(nd);
225 CheckLabels2(milab);
226 c[*n].SetLabels(milab);
a86176e3 227
228 if(repa->GetUseUnfoldingInClusterFinderSSD()==kTRUE) {
d695268b 229 // Note: fUseUnfoldingInClusterFinderSSD=kFALSE by default in RecoParam
230
a86176e3 231 //Split suspiciously big cluster
232 if (nd>4&&nd<25) {
233 c[*n].SetY(y/q-0.25*nd);
234 c[*n].SetQ(0.5*q);
235 (*n)++;
236 if (*n==kMax) {
237 Error("FindClustersSSD","Too many 1D clusters !");
238 return;
239 }
240 c[*n].SetY(y/q+0.25*nd);
241 c[*n].SetQ(0.5*q);
242 c[*n].SetNd(nd);
243 c[*n].SetLabels(milab);
244 }
245
246 } // unfolding is on
247
04366a57 248 (*n)++;
249 if (*n==kMax) {
250 Error("FindClustersSSD","Too many 1D clusters !");
251 return;
252 }
d695268b 253
254 } // flag5 set
255
256 // reset everything
04366a57 257 y=q=qmax=0.;
258 nd=0;
d695268b 259 flag5=0;
04366a57 260 lab[0]=lab[1]=lab[2]=-2;
d695268b 261 for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2;
262
263 // if side changed from P to N, switch to pos 1D clusters
264 // (if for some reason the side changed from N to P then do the opposite)
265 if (flag!=d->GetCoord1())
266 { if(!flag) {n=&np; c=pos;} else {n=&nn; c=neg;} }
267
268 } // end create new 1D cluster from previous neighboring digits
269
270 // continues adding digits to the previous cluster
271 // or start a new one
04366a57 272 flag=d->GetCoord1();
273 q += d->GetSignal();
274 y += d->GetCoord2()*d->GetSignal();
275 nd++;
d695268b 276
277 if(d->IsSideP()) {
278 noise = cal->GetNoiseP(d->GetStripNumber());
279 gain = cal->GetGainP(d->GetStripNumber());
280 }
281 else {
282 noise = cal->GetNoiseN(d->GetStripNumber());
283 gain = cal->GetGainN(d->GetStripNumber());
284 }
285 noise*=gain;
286 noise=cal->ADCToKeV(noise); // converts the charge in KeV from ADC units
287
288 if(d->GetSignal()>fgkThreshold*noise) flag5=1;
289
290 /*
291 cout<<d->GetSignal()<<" "<<noise<<" "<<flag5<<" "<<
292 d->GetCoord1()<<" "<<d->GetCoord2()<<endl;
293 */
294
04366a57 295 if (d->GetSignal()>qmax) {
296 qmax=d->GetSignal();
297 lab[0]=d->GetTrack(0); lab[1]=d->GetTrack(1); lab[2]=d->GetTrack(2);
298 }
299 for (Int_t ilab=0;ilab<10;ilab++) {
300 if (d->GetTrack(ilab)>=0) AddLabel(milab, (d->GetTrack(ilab)));
301 }
302 curr=strip;
a86176e3 303
d695268b 304
305 } // loop over digits, no more digits in the digits array
306
307
308 // add the last 1D cluster
309 if(flag5) {
310
311 // cout<<"here2"<<endl;
312
313 c[*n].SetY(y/q);
314 c[*n].SetQ(q);
315 c[*n].SetNd(nd);
316 c[*n].SetLabels(lab);
a86176e3 317
d695268b 318 if(repa->GetUseUnfoldingInClusterFinderSSD()==kTRUE) {
319
320 //Split suspiciously big cluster
321 if (nd>4 && nd<25) {
322 c[*n].SetY(y/q-0.25*nd);
323 c[*n].SetQ(0.5*q);
324 (*n)++;
325 if (*n==kMax) {
326 Error("FindClustersSSD","Too many 1D clusters !");
327 return;
328 }
329 c[*n].SetY(y/q+0.25*nd);
330 c[*n].SetQ(0.5*q);
331 c[*n].SetNd(nd);
332 c[*n].SetLabels(lab);
a86176e3 333 }
d695268b 334 } // unfolding is on
335
336 (*n)++;
337 if (*n==kMax) {
338 Error("FindClustersSSD","Too many 1D clusters !");
339 return;
a86176e3 340 }
04366a57 341
d695268b 342 } // if flag5 last 1D cluster added
343
344
345 //------------------------------------------------------
346 // call FindClustersSSD to pair neg and pos 1D clusters
347 // and create recpoints from the crosses
348 // Note1: neg are Pside and pos are Nside!!
349 // Note2: if there are no Pside digits nn=0 (bad strips!!) (same for Nside)
350 //
351 // cout<<nn<<" Pside and "<<np<<" Nside clusters"<<endl;
04366a57 352 FindClustersSSD(neg, nn, pos, np);
d695268b 353 //
354 //-----------------------------------------------------
355
04366a57 356}
357
358
359void AliITSClusterFinderV2SSD::RawdataToClusters(AliRawReader* rawReader,TClonesArray** clusters){
360
361 //------------------------------------------------------------
362 // This function creates ITS clusters from raw data
363 //------------------------------------------------------------
364 rawReader->Reset();
365 AliITSRawStreamSSD inputSSD(rawReader);
366 FindClustersSSD(&inputSSD,clusters);
367
368}
369
3a4139a2 370void AliITSClusterFinderV2SSD::FindClustersSSD(AliITSRawStreamSSD* input,
04366a57 371 TClonesArray** clusters)
372{
373 //------------------------------------------------------------
374 // Actual SSD cluster finder for raw data
375 //------------------------------------------------------------
a86176e3 376
377 static AliITSRecoParam *repa = NULL;
378 if(!repa){
379 repa = (AliITSRecoParam*) AliITSReconstructor::GetRecoParam();
380 if(!repa){
ed446fa3 381 repa = AliITSRecoParam::GetHighFluxParam();
a86176e3 382 AliWarning("Using default AliITSRecoParam class");
383 }
384 }
385
04366a57 386 Int_t nClustersSSD = 0;
387 const Int_t kMax = 1000;
388 Ali1Dcluster clusters1D[2][kMax];
389 Int_t nClusters[2] = {0, 0};
390 Int_t lab[3]={-2,-2,-2};
391 Float_t q = 0.;
392 Float_t y = 0.;
393 Int_t nDigits = 0;
3a4139a2 394 Float_t gain=0;
308b5ea4 395 Float_t noise=0.;
5ffe058f 396 // Float_t pedestal=0.;
308b5ea4 397 Float_t oldnoise=0.;
b4704be3 398 AliITSCalibrationSSD* cal=NULL;
308b5ea4 399
308b5ea4 400 Int_t matrix[12][1536];
d5ad0697 401 Int_t iddl=-1;
402 Int_t iad=-1;
403 Int_t oddl = -1;
404 Int_t oad = -1;
405 Int_t oadc = -1;
406 Int_t ostrip = -1;
308b5ea4 407 Int_t osignal = 65535;
308b5ea4 408 Int_t n=0;
409 Bool_t next=0;
04366a57 410
411 // read raw data input stream
412 while (kTRUE) {
308b5ea4 413
414 // reset signal matrix
415 for(Int_t i=0; i<12; i++) { for(Int_t j=0; j<1536; j++) { matrix[i][j] = 65535;} }
b42cfa25 416
cbf26d57 417 if((osignal!=65535)&&(ostrip>0)&&(ostrip<1536)) {
308b5ea4 418 n++;
419 matrix[oadc][ostrip] = osignal; // recover data from previous occurence of input->Next()
420 }
b42cfa25 421
308b5ea4 422 // buffer data for ddl=iddl and ad=iad
423 while(kTRUE) {
bc4dd89a 424
308b5ea4 425 next = input->Next();
d5ad0697 426 if((!next)&&(input->flag)) continue;
427 Int_t ddl=input->GetDDL();
428 Int_t ad=input->GetAD();
429 Int_t adc = input->GetADC(); adc = (adc<6)? adc : adc - 2;
430 Int_t strip = input->GetStrip();
308b5ea4 431 if(input->GetSideFlag()) strip=1535-strip;
d5ad0697 432 Int_t signal = input->GetSignal();
d036ccd3 433
cbf26d57 434 if((ddl==iddl)&&(ad==iad)&&(strip>0)&&(strip<1536)) {n++; matrix[adc][strip] = signal;}
4ae225e9 435 else {if ((strip<1536) && (strip>0)) {oddl=iddl; oad=iad; oadc = adc; ostrip = strip; osignal=signal; iddl=ddl; iad=ad; break;}}
308b5ea4 436
b42cfa25 437 if(!next) {oddl=iddl; oad=iad; oadc = adc; ostrip = strip; osignal=signal; iddl=ddl; iad=ad; break;}
438 //break;
308b5ea4 439 }
b42cfa25 440
441 // No SSD data
d5ad0697 442 if(!next && oddl<0) break;
b42cfa25 443
308b5ea4 444 if(n==0) continue; // first occurence
d036ccd3 445 n=0; //osignal=0;
42ed6062 446
447 Float_t dStrip = 0;
448 if (repa->GetUseCosmicRunShiftsSSD()) { // Special condition for 2007/2008 cosmic data
449 dStrip = fgkCosmic2008StripShifts[oddl][oad-1];
450 }
b6b7d649 451 if (TMath::Abs(dStrip) > 1.5)
452 AliError(Form("Indexing error ? oddl = %d, dStrip %f\n",oddl,dStrip));
308b5ea4 453 // fill 1Dclusters
454 for(Int_t iadc=0; iadc<12; iadc++) { // loop over ADC index for ddl=oddl and ad=oad
b42cfa25 455
308b5ea4 456 Int_t iimod = (oad - 1) * 12 + iadc;
d5ad0697 457 Int_t iModule = AliITSRawStreamSSD::GetModuleNumber(oddl,iimod);
308b5ea4 458 if(iModule==-1) continue;
308b5ea4 459 cal = (AliITSCalibrationSSD*)GetResp(iModule);
42ed6062 460
308b5ea4 461 Bool_t first = 0;
d695268b 462 Bool_t flag5 = 0;
b42cfa25 463
bc4dd89a 464 /*
308b5ea4 465 for(Int_t istrip=0; istrip<768; istrip++) { // P-side
d5ad0697 466 Int_t signal = matrix[iadc][istrip];
b42cfa25 467 pedestal = cal->GetPedestalP(istrip);
e7f0e76b 468 matrix[iadc][istrip]=signal-(Int_t)pedestal;
bc4dd89a 469 }
470 */
471
472 /*
b42cfa25 473 Float_t cmode=0;
474 for(Int_t l=0; l<6; l++) {
475 cmode=0;
476 for(Int_t n=20; n<108; n++) cmode+=matrix[iadc][l*128+n];
477 cmode/=88.;
e7f0e76b 478 for(Int_t n=0; n<128; n++) matrix[iadc][l*128+n]-=(Int_t)cmode;
b42cfa25 479
480 }
bc4dd89a 481 */
482
a64f9843 483 Int_t istrip=0;
484 for(istrip=0; istrip<768; istrip++) { // P-side
b42cfa25 485
486 Int_t signal = TMath::Abs(matrix[iadc][istrip]);
487
308b5ea4 488 oldnoise = noise;
b42cfa25 489 noise = cal->GetNoiseP(istrip); if(noise<1.) signal = 65535;
a64f9843 490 if(signal<3*noise) signal = 65535; // in case ZS was not done in hw do it now
491
ced4d9bc 492 // if(cal->IsPChannelBad(istrip)) signal=0;
308b5ea4 493
494 if (signal!=65535) {
495 gain = cal->GetGainP(istrip);
496 signal = (Int_t) ( signal * gain ); // signal is corrected for gain
d695268b 497 if(signal>fgkThreshold*noise) flag5=1;
308b5ea4 498 signal = (Int_t) cal->ADCToKeV( signal ); // signal is converted in KeV
b42cfa25 499
308b5ea4 500 q += signal; // add digit to current cluster
501 y += istrip * signal;
502 nDigits++;
503 first=1;
504 }
b42cfa25 505
308b5ea4 506 else if(first) {
b42cfa25 507
d695268b 508 if ( (nDigits>0) && flag5 ) {
b42cfa25 509
308b5ea4 510 Ali1Dcluster& cluster = clusters1D[0][nClusters[0]++];
a64f9843 511
42ed6062 512 if(q!=0) cluster.SetY(y/q + dStrip);
513 else cluster.SetY(istrip + dStrip -1);
a64f9843 514
308b5ea4 515 cluster.SetQ(q);
516 cluster.SetNd(nDigits);
517 cluster.SetLabels(lab);
518
a86176e3 519 if(repa->GetUseUnfoldingInClusterFinderSSD()==kTRUE) {
520
521 //Split suspiciously big cluster
522 if (nDigits > 4&&nDigits < 25) {
42ed6062 523 if(q!=0) cluster.SetY(y/q + dStrip - 0.25*nDigits);
524 else cluster.SetY(istrip-1 + dStrip - 0.25*nDigits);
a86176e3 525 cluster.SetQ(0.5*q);
526 if (nClusters[0] == kMax) {
527 Error("FindClustersSSD", "Too many 1D clusters !");
528 return;
529 }
530 Ali1Dcluster& cluster2 = clusters1D[0][nClusters[0]++];
42ed6062 531 if(q!=0) cluster2.SetY(y/q + dStrip + 0.25*nDigits);
532 else cluster2.SetY(istrip-1 + dStrip + 0.25*nDigits);
a86176e3 533 cluster2.SetQ(0.5*q);
534 cluster2.SetNd(nDigits);
535 cluster2.SetLabels(lab);
308b5ea4 536 }
a86176e3 537 } // unfolding is on
308b5ea4 538 }
539
540 y = q = 0.;
541 nDigits = 0;
542 first=0;
d695268b 543 flag5=0;
04366a57 544 }
a86176e3 545
308b5ea4 546 } // loop over strip on P-side
a86176e3 547
308b5ea4 548 // if last strip does have signal
549 if(first) {
550
d695268b 551 if ( (nDigits>0) && flag5 ) {
a86176e3 552
553 Ali1Dcluster& cluster = clusters1D[0][nClusters[0]++];
a64f9843 554
42ed6062 555 if(q!=0) cluster.SetY(y/q + dStrip);
556 else cluster.SetY(istrip - 1 + dStrip);
a64f9843 557
a86176e3 558 cluster.SetQ(q);
559 cluster.SetNd(nDigits);
560 cluster.SetLabels(lab);
561
562 if(repa->GetUseUnfoldingInClusterFinderSSD()==kTRUE) {
308b5ea4 563
564 //Split suspiciously big cluster
565 if (nDigits > 4&&nDigits < 25) {
42ed6062 566 if(q!=0) cluster.SetY(y/q + dStrip - 0.25*nDigits);
567 else cluster.SetY(istrip-1 + dStrip - 0.25*nDigits);
308b5ea4 568 cluster.SetQ(0.5*q);
569 if (nClusters[0] == kMax) {
570 Error("FindClustersSSD", "Too many 1D clusters !");
571 return;
572 }
573 Ali1Dcluster& cluster2 = clusters1D[0][nClusters[0]++];
42ed6062 574 if(q!=0) cluster2.SetY(y/q + dStrip + 0.25*nDigits);
575 else cluster2.SetY(istrip-1 + dStrip + 0.25*nDigits);
308b5ea4 576 cluster2.SetQ(0.5*q);
577 cluster2.SetNd(nDigits);
578 cluster2.SetLabels(lab);
579 }
a86176e3 580 } // unfolding is on
581
582 }
583 y = q = 0.;
584 nDigits = 0;
585 first=0;
d695268b 586 flag5=0;
04366a57 587 }
308b5ea4 588
bc4dd89a 589 /*
b42cfa25 590 for(Int_t istrip=768; istrip<1536; istrip++) { // P-side
591 Int_t signal = matrix[iadc][istrip];
592 pedestal = cal->GetPedestalN(1535-istrip);
e7f0e76b 593 matrix[iadc][istrip]=signal-(Int_t)pedestal;
b42cfa25 594 }
bc4dd89a 595 */
b42cfa25 596
bc4dd89a 597 /*
b42cfa25 598 for(Int_t l=6; l<12; l++) {
599 Float_t cmode=0;
600 for(Int_t n=20; n<108; n++) cmode+=matrix[iadc][l*128+n];
601 cmode/=88.;
e7f0e76b 602 for(Int_t n=0; n<128; n++) matrix[iadc][l*128+n]-=(Int_t)cmode;
b42cfa25 603 }
bc4dd89a 604 */
b42cfa25 605
308b5ea4 606 oldnoise = 0.;
607 noise = 0.;
a64f9843 608 Int_t strip=0;
4adcf390 609 for(Int_t iistrip=768; iistrip<1536; iistrip++) { // N-side
308b5ea4 610
4adcf390 611 Int_t signal = TMath::Abs(matrix[iadc][iistrip]);
4adcf390 612 strip = 1535-iistrip;
b42cfa25 613
308b5ea4 614 oldnoise = noise;
b42cfa25 615 noise = cal->GetNoiseN(strip); if(noise<1.) signal=65535;
a64f9843 616
d036ccd3 617 // if(cal->IsNChannelBad(strip)) signal=0;
a64f9843 618
619 if(signal<3*noise) signal = 65535; // in case ZS was not done in hw do it now
308b5ea4 620
621 if (signal!=65535) {
622 gain = cal->GetGainN(strip);
623 signal = (Int_t) ( signal * gain); // signal is corrected for gain
d695268b 624 if(signal>fgkThreshold*noise) flag5=1;
308b5ea4 625 signal = (Int_t) cal->ADCToKeV( signal ); // signal is converted in KeV
626
627 // add digit to current cluster
628 q += signal;
629 y += strip * signal;
630 nDigits++;
631 first=1;
632 }
633
634 else if(first) {
a86176e3 635
d695268b 636 if ( (nDigits>0) && flag5 ) {
308b5ea4 637
638 Ali1Dcluster& cluster = clusters1D[1][nClusters[1]++];
a64f9843 639
42ed6062 640 if(q!=0) cluster.SetY(y/q - dStrip);
641 else cluster.SetY(strip+1 - dStrip);
a64f9843 642
308b5ea4 643 cluster.SetQ(q);
644 cluster.SetNd(nDigits);
645 cluster.SetLabels(lab);
646
a86176e3 647 if(repa->GetUseUnfoldingInClusterFinderSSD()==kTRUE) {
648
649 //Split suspiciously big cluster
650 if (nDigits > 4&&nDigits < 25) {
42ed6062 651 cluster.SetY(y/q - dStrip - 0.25*nDigits);
a86176e3 652 cluster.SetQ(0.5*q);
653 if (nClusters[1] == kMax) {
654 Error("FindClustersSSD", "Too many 1D clusters !");
655 return;
656 }
657 Ali1Dcluster& cluster2 = clusters1D[1][nClusters[1]++];
42ed6062 658 cluster2.SetY(y/q - dStrip + 0.25*nDigits);
a86176e3 659 cluster2.SetQ(0.5*q);
660 cluster2.SetNd(nDigits);
661 cluster2.SetLabels(lab);
662 }
663 } // unfolding is on
664 }
665
308b5ea4 666 y = q = 0.;
667 nDigits = 0;
668 first=0;
d695268b 669 flag5=0;
308b5ea4 670 }
671
672 } // loop over strips on N-side
04366a57 673
308b5ea4 674 if(first) {
675
d695268b 676 if ( (nDigits>0) && flag5 ) {
a86176e3 677
678 Ali1Dcluster& cluster = clusters1D[1][nClusters[1]++];
a64f9843 679
42ed6062 680 if(q!=0) cluster.SetY(y/q - dStrip);
681 else cluster.SetY(strip - dStrip + 1);
a64f9843 682
a86176e3 683 cluster.SetQ(q);
684 cluster.SetNd(nDigits);
685 cluster.SetLabels(lab);
686
687 if(repa->GetUseUnfoldingInClusterFinderSSD()==kTRUE) {
308b5ea4 688
689 //Split suspiciously big cluster
690 if (nDigits > 4&&nDigits < 25) {
42ed6062 691 if(q!=0) cluster.SetY(y/q - dStrip - 0.25*nDigits);
692 else cluster.SetY(strip+1 - dStrip - 0.25*nDigits);
308b5ea4 693 cluster.SetQ(0.5*q);
694 if (nClusters[1] == kMax) {
695 Error("FindClustersSSD", "Too many 1D clusters !");
696 return;
697 }
698 Ali1Dcluster& cluster2 = clusters1D[1][nClusters[1]++];
42ed6062 699 if(q!=0) cluster2.SetY(y/q - dStrip + 0.25*nDigits);
700 else cluster2.SetY(strip+1 - dStrip + 0.25*nDigits);
308b5ea4 701 cluster2.SetQ(0.5*q);
702 cluster2.SetNd(nDigits);
703 cluster2.SetLabels(lab);
704 }
a86176e3 705 } // unfolding is on
706 }
707
708 y = q = 0.;
709 nDigits = 0;
710 first=0;
d695268b 711 flag5=0;
308b5ea4 712 }
713
714 // create recpoints
715 if((nClusters[0])&&(nClusters[1])) {
a86176e3 716
00a7cc50 717 clusters[iModule] = new TClonesArray("AliITSRecPoint");
04366a57 718 fModule = iModule;
719 FindClustersSSD(&clusters1D[0][0], nClusters[0],
720 &clusters1D[1][0], nClusters[1], clusters[iModule]);
4adcf390 721 Int_t nClustersn = clusters[iModule]->GetEntriesFast();
722 nClustersSSD += nClustersn;
04366a57 723 }
724
04366a57 725 nClusters[0] = nClusters[1] = 0;
726 y = q = 0.;
727 nDigits = 0;
3a4139a2 728
308b5ea4 729 } // loop over adc
04366a57 730
d5ad0697 731 if(!next) break;
04366a57 732 }
308b5ea4 733
04366a57 734 Info("FindClustersSSD", "found clusters in ITS SSD: %d", nClustersSSD);
735}
736
737void AliITSClusterFinderV2SSD::
738FindClustersSSD(Ali1Dcluster* neg, Int_t nn,
739 Ali1Dcluster* pos, Int_t np,
740 TClonesArray *clusters) {
741 //------------------------------------------------------------
742 // Actual SSD cluster finder
743 //------------------------------------------------------------
b4704be3 744
745 const TGeoHMatrix *mT2L=AliITSgeomTGeo::GetTracking2LocalMatrix(fModule);
746
7101948c 747 //---------------------------------------
748 // load recoparam
749 //
750 static AliITSRecoParam *repa = NULL;
751 if(!repa){
752 repa = (AliITSRecoParam*) AliITSReconstructor::GetRecoParam();
753 if(!repa){
754 repa = AliITSRecoParam::GetHighFluxParam();
755 AliWarning("Using default AliITSRecoParam class");
756 }
757 }
758
04366a57 759 TClonesArray &cl=*clusters;
d036ccd3 760
8be4e1b1 761 AliITSsegmentationSSD *seg = dynamic_cast<AliITSsegmentationSSD*>(fDetTypeRec->GetSegmentationModel(2));
762 if (fModule>fLastSSD1)
763 seg->SetLayer(6);
764 else
765 seg->SetLayer(5);
766
767 Float_t hwSSD = seg->Dx()*1e-4/2;
768 Float_t hlSSD = seg->Dz()*1e-4/2;
769
04366a57 770 Int_t idet=fNdet[fModule];
771 Int_t ncl=0;
d036ccd3 772
04366a57 773 //
8be4e1b1 774 Int_t *cnegative = new Int_t[np];
775 Int_t *cused1 = new Int_t[np];
776 Int_t *negativepair = new Int_t[10*np];
777 Int_t *cpositive = new Int_t[nn];
778 Int_t *cused2 = new Int_t[nn];
779 Int_t *positivepair = new Int_t[10*nn];
780 for (Int_t i=0;i<np;i++) {cnegative[i]=0; cused1[i]=0;}
781 for (Int_t i=0;i<nn;i++) {cpositive[i]=0; cused2[i]=0;}
782 for (Int_t i=0;i<10*np;i++) {negativepair[i]=0;}
783 for (Int_t i=0;i<10*nn;i++) {positivepair[i]=0;}
c157c94e 784
308b5ea4 785 if ((np*nn) > fgPairsSize) {
d036ccd3 786
308b5ea4 787 if (fgPairs) delete [] fgPairs;
788 fgPairsSize = 4*np*nn;
789 fgPairs = new Short_t[fgPairsSize];
c157c94e 790 }
308b5ea4 791 memset(fgPairs,0,sizeof(Short_t)*np*nn);
792
04366a57 793 //
794 // find available pairs
795 //
796 for (Int_t i=0; i<np; i++) {
d036ccd3 797 Float_t yp=pos[i].GetY();
a64f9843 798 if ( (pos[i].GetQ()>0) && (pos[i].GetQ()<3) ) continue;
04366a57 799 for (Int_t j=0; j<nn; j++) {
a64f9843 800 if ( (neg[j].GetQ()>0) && (neg[j].GetQ()<3) ) continue;
d036ccd3 801 Float_t yn=neg[j].GetY();
802
8be4e1b1 803 Float_t xt, zt;
804 seg->GetPadCxz(yn, yp, xt, zt);
d695268b 805 //cout<<yn<<" "<<yp<<" "<<xt<<" "<<zt<<endl;
d036ccd3 806
8be4e1b1 807 if (TMath::Abs(xt)<hwSSD+0.01)
808 if (TMath::Abs(zt)<hlSSD+0.01*(neg[j].GetNd()+pos[i].GetNd())) {
710f576f 809 Int_t in = i*10+cnegative[i];
810 Int_t ip = j*10+cpositive[j];
811 if ((in < 10*np) && (ip < 10*nn)) {
812 negativepair[in] =j; //index
813 positivepair[ip] =i;
814 cnegative[i]++; //counters
815 cpositive[j]++;
816 fgPairs[i*nn+j]=100;
817 }
818 else
819 AliError(Form("Index out of range: ip=%d, in=%d",ip,in));
04366a57 820 }
821 }
822 }
308b5ea4 823
824 //
825 // try to recover points out of but close to the module boundaries
04366a57 826 //
827 for (Int_t i=0; i<np; i++) {
d036ccd3 828 Float_t yp=pos[i].GetY();
a64f9843 829 if ( (pos[i].GetQ()>0) && (pos[i].GetQ()<3) ) continue;
04366a57 830 for (Int_t j=0; j<nn; j++) {
a64f9843 831 if ( (neg[j].GetQ()>0) && (neg[j].GetQ()<3) ) continue;
308b5ea4 832 // if both 1Dclusters have an other cross continue
04366a57 833 if (cpositive[j]&&cnegative[i]) continue;
8be4e1b1 834 Float_t yn=neg[j].GetY();
d036ccd3 835
8be4e1b1 836 Float_t xt, zt;
837 seg->GetPadCxz(yn, yp, xt, zt);
d036ccd3 838
8be4e1b1 839 if (TMath::Abs(xt)<hwSSD+0.1)
840 if (TMath::Abs(zt)<hlSSD+0.15) {
308b5ea4 841 // tag 1Dcluster (eventually will produce low quality recpoint)
04366a57 842 if (cnegative[i]==0) pos[i].SetNd(100); // not available pair
843 if (cpositive[j]==0) neg[j].SetNd(100); // not available pair
710f576f 844 Int_t in = i*10+cnegative[i];
845 Int_t ip = j*10+cpositive[j];
846 if ((in < 10*np) && (ip < 10*nn)) {
847 negativepair[in] =j; //index
848 positivepair[ip] =i;
849 cnegative[i]++; //counters
850 cpositive[j]++;
851 fgPairs[i*nn+j]=100;
852 }
853 else
854 AliError(Form("Index out of range: ip=%d, in=%d",ip,in));
04366a57 855 }
856 }
857 }
308b5ea4 858
04366a57 859 //
d695268b 860 Float_t lp[6];
04366a57 861 Int_t milab[10];
862 Double_t ratio;
863
b42cfa25 864
a86176e3 865 if(repa->GetUseChargeMatchingInClusterFinderSSD()==kTRUE) {
308b5ea4 866
867
a86176e3 868 //
869 // sign gold tracks
870 //
871 for (Int_t ip=0;ip<np;ip++){
8be4e1b1 872 Float_t xbest=1000,zbest=1000,qbest=0;
a86176e3 873 //
874 // select gold clusters
875 if ( (cnegative[ip]==1) && cpositive[negativepair[10*ip]]==1){
d036ccd3 876 Float_t yp=pos[ip].GetY();
a86176e3 877 Int_t j = negativepair[10*ip];
a64f9843 878
879 if( (pos[ip].GetQ()==0) && (neg[j].GetQ() ==0) ) {
880 // both bad, hence continue;
881 // mark both as used (to avoid recover at the end)
882 cused1[ip]++;
883 cused2[j]++;
884 continue;
885 }
886
a86176e3 887 ratio = (pos[ip].GetQ()-neg[j].GetQ())/(pos[ip].GetQ()+neg[j].GetQ());
d695268b 888 //cout<<"ratio="<<ratio<<endl;
a64f9843 889
890 // charge matching (note that if posQ or negQ is 0 -> ratio=1 and the following condition is met
2069484c 891 if (TMath::Abs(ratio)>0.2) continue; // note: 0.2=3xsigma_ratio calculated in cosmics tests
a64f9843 892
a86176e3 893 //
8be4e1b1 894 Float_t yn=neg[j].GetY();
895
896 Float_t xt, zt;
897 seg->GetPadCxz(yn, yp, xt, zt);
898
899 xbest=xt; zbest=zt;
2069484c 900
2069484c 901
a86176e3 902 qbest=0.5*(pos[ip].GetQ()+neg[j].GetQ());
d695268b 903 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 904
a86176e3 905 {
8be4e1b1 906 Double_t loc[3]={xbest,0.,zbest},trk[3]={0.,0.,0.};
a86176e3 907 mT2L->MasterToLocal(loc,trk);
908 lp[0]=trk[1];
909 lp[1]=trk[2];
00a7cc50 910 }
a86176e3 911 lp[4]=qbest; //Q
912 for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2;
913 for (Int_t ilab=0;ilab<3;ilab++){
914 milab[ilab] = pos[ip].GetLabel(ilab);
915 milab[ilab+3] = neg[j].GetLabel(ilab);
00a7cc50 916 }
04366a57 917 //
a86176e3 918 CheckLabels2(milab);
919 milab[3]=(((ip<<10) + j)<<10) + idet; // pos|neg|det
920 Int_t info[3] = {pos[ip].GetNd(),neg[j].GetNd(),fNlayer[fModule]};
d695268b 921
922 lp[2]=0.0022*0.0022; //SigmaY2
923 lp[3]=0.110*0.110; //SigmaZ2
924 // out-of-diagonal element of covariance matrix
925 if( (info[0]==1) && (info[1]==1) ) lp[5]=-0.00012;
926 else if ( (info[0]>1) && (info[1]>1) ) {
927 lp[2]=0.0016*0.0016; //SigmaY2
928 lp[3]=0.08*0.08; //SigmaZ2
929 lp[5]=-0.00006;
930 }
931 else {
7101948c 932 lp[3]=0.093*0.093;
d695268b 933 if (info[0]==1) { lp[5]=-0.00014;}
934 else { lp[2]=0.0017*0.0017; lp[5]=-0.00004;}
935 }
936
a86176e3 937 AliITSRecPoint * cl2;
938
939 if(clusters){ // Note clusters != 0 when method is called for rawdata
04366a57 940
a86176e3 941
942 cl2 = new (cl[ncl]) AliITSRecPoint(milab,lp,info);
943
a86176e3 944 cl2->SetChargeRatio(ratio);
945 cl2->SetType(1);
946 fgPairs[ip*nn+j]=1;
a64f9843 947
a86176e3 948 if ((pos[ip].GetNd()+neg[j].GetNd())>6){ //multi cluster
949 cl2->SetType(2);
950 fgPairs[ip*nn+j]=2;
04366a57 951 }
a64f9843 952
953 if(pos[ip].GetQ()==0) cl2->SetType(3);
7b208ef4 954 if(neg[j].GetQ()==0) cl2->SetType(4);
a64f9843 955
a86176e3 956 cused1[ip]++;
957 cused2[j]++;
958
959 }
960 else{ // Note clusters == 0 when method is called for digits
961
962 cl2 = new AliITSRecPoint(milab,lp,info);
963
a86176e3 964 cl2->SetChargeRatio(ratio);
965 cl2->SetType(1);
966 fgPairs[ip*nn+j]=1;
a64f9843 967
a86176e3 968 if ((pos[ip].GetNd()+neg[j].GetNd())>6){ //multi cluster
969 cl2->SetType(2);
970 fgPairs[ip*nn+j]=2;
00a7cc50 971 }
a64f9843 972
973 if(pos[ip].GetQ()==0) cl2->SetType(3);
7b208ef4 974 if(neg[j].GetQ()==0) cl2->SetType(4);
a64f9843 975
a86176e3 976 cused1[ip]++;
977 cused2[j]++;
a64f9843 978
a86176e3 979 fDetTypeRec->AddRecPoint(*cl2);
980 }
981 ncl++;
982 }
983 }
984
985 for (Int_t ip=0;ip<np;ip++){
8be4e1b1 986 Float_t xbest=1000,zbest=1000,qbest=0;
a86176e3 987 //
988 //
989 // select "silber" cluster
990 if ( cnegative[ip]==1 && cpositive[negativepair[10*ip]]==2){
991 Int_t in = negativepair[10*ip];
992 Int_t ip2 = positivepair[10*in];
993 if (ip2==ip) ip2 = positivepair[10*in+1];
994 Float_t pcharge = pos[ip].GetQ()+pos[ip2].GetQ();
a64f9843 995
d695268b 996
997
998 ratio = (pcharge-neg[in].GetQ())/(pcharge+neg[in].GetQ());
999 if ( (TMath::Abs(ratio)<0.2) && (pcharge!=0) ) {
1000 //if ( (TMath::Abs(pcharge-neg[in].GetQ())<30) && (pcharge!=0) ) { //
a64f9843 1001
a86176e3 1002 //
1003 // add first pair
a64f9843 1004 if ( (fgPairs[ip*nn+in]==100)&&(pos[ip].GetQ() ) ) { //
1005
d036ccd3 1006 Float_t yp=pos[ip].GetY();
8be4e1b1 1007 Float_t yn=neg[in].GetY();
1008
1009 Float_t xt, zt;
1010 seg->GetPadCxz(yn, yp, xt, zt);
1011
1012 xbest=xt; zbest=zt;
2069484c 1013
a86176e3 1014 qbest =pos[ip].GetQ();
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];
1019
a86176e3 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[in].GetLabel(ilab);
1025 }
1026 //
1027 CheckLabels2(milab);
1028 ratio = (pos[ip].GetQ()-neg[in].GetQ())/(pos[ip].GetQ()+neg[in].GetQ());
1029 milab[3]=(((ip<<10) + in)<<10) + idet; // pos|neg|det
1030 Int_t info[3] = {pos[ip].GetNd(),neg[in].GetNd(),fNlayer[fModule]};
1031
d695268b 1032 lp[2]=0.0022*0.0022; //SigmaY2
1033 lp[3]=0.110*0.110; //SigmaZ2
1034 // out-of-diagonal element of covariance matrix
1035 if( (info[0]==1) && (info[1]==1) ) lp[5]=-0.00012;
1036 else if ( (info[0]>1) && (info[1]>1) ) {
1037 lp[2]=0.0016*0.0016; //SigmaY2
1038 lp[3]=0.08*0.08; //SigmaZ2
1039 lp[5]=-0.00006;
1040 }
1041 else {
7101948c 1042 lp[3]=0.093*0.093;
d695268b 1043 if (info[0]==1) { lp[5]=-0.00014;}
1044 else { lp[2]=0.0017*0.0017; lp[5]=-0.00004;}
1045 }
1046
a86176e3 1047 AliITSRecPoint * cl2;
1048 if(clusters){
1049
1050 cl2 = new (cl[ncl]) AliITSRecPoint(milab,lp,info);
a86176e3 1051 cl2->SetChargeRatio(ratio);
1052 cl2->SetType(5);
1053 fgPairs[ip*nn+in] = 5;
1054 if ((pos[ip].GetNd()+neg[in].GetNd())>6){ //multi cluster
1055 cl2->SetType(6);
1056 fgPairs[ip*nn+in] = 6;
1057 }
1058 }
1059 else{
1060 cl2 = new AliITSRecPoint(milab,lp,info);
1061 cl2->SetChargeRatio(ratio);
1062 cl2->SetType(5);
1063 fgPairs[ip*nn+in] = 5;
1064 if ((pos[ip].GetNd()+neg[in].GetNd())>6){ //multi cluster
1065 cl2->SetType(6);
1066 fgPairs[ip*nn+in] = 6;
1067 }
a86176e3 1068
1069 fDetTypeRec->AddRecPoint(*cl2);
1070 }
1071 ncl++;
04366a57 1072 }
04366a57 1073
a64f9843 1074
04366a57 1075 //
a86176e3 1076 // add second pair
00a7cc50 1077
a86176e3 1078 // if (!(cused1[ip2] || cused2[in])){ //
a64f9843 1079 if ( (fgPairs[ip2*nn+in]==100) && (pos[ip2].GetQ()) ) {
1080
d036ccd3 1081 Float_t yp=pos[ip2].GetY();
8be4e1b1 1082 Float_t yn=neg[in].GetY();
1083
1084 Float_t xt, zt;
1085 seg->GetPadCxz(yn, yp, xt, zt);
1086
1087 xbest=xt; zbest=zt;
2069484c 1088
a86176e3 1089 qbest =pos[ip2].GetQ();
a64f9843 1090
8be4e1b1 1091 Double_t loc[3]={xbest,0.,zbest},trk[3]={0.,0.,0.};
a64f9843 1092 mT2L->MasterToLocal(loc,trk);
1093 lp[0]=trk[1];
1094 lp[1]=trk[2];
1095
a86176e3 1096 lp[4]=qbest; //Q
1097 for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2;
1098 for (Int_t ilab=0;ilab<3;ilab++){
1099 milab[ilab] = pos[ip2].GetLabel(ilab);
1100 milab[ilab+3] = neg[in].GetLabel(ilab);
00a7cc50 1101 }
a86176e3 1102 //
1103 CheckLabels2(milab);
1104 ratio = (pos[ip2].GetQ()-neg[in].GetQ())/(pos[ip2].GetQ()+neg[in].GetQ());
1105 milab[3]=(((ip2<<10) + in)<<10) + idet; // pos|neg|det
1106 Int_t info[3] = {pos[ip2].GetNd(),neg[in].GetNd(),fNlayer[fModule]};
308b5ea4 1107
d695268b 1108 lp[2]=0.0022*0.0022; //SigmaY2
1109 lp[3]=0.110*0.110; //SigmaZ2
1110 // out-of-diagonal element of covariance matrix
1111 if( (info[0]==1) && (info[1]==1) ) lp[5]=-0.00012;
1112 else if ( (info[0]>1) && (info[1]>1) ) {
1113 lp[2]=0.0016*0.0016; //SigmaY2
1114 lp[3]=0.08*0.08; //SigmaZ2
1115 lp[5]=-0.00006;
1116 }
1117 else {
7101948c 1118 lp[3]=0.093*0.093;
d695268b 1119 if (info[0]==1) { lp[5]=-0.00014;}
1120 else { lp[2]=0.0017*0.0017; lp[5]=-0.00004;}
1121 }
1122
a86176e3 1123 AliITSRecPoint * cl2;
1124 if(clusters){
1125 cl2 = new (cl[ncl]) AliITSRecPoint(milab,lp,info);
1126
a86176e3 1127 cl2->SetChargeRatio(ratio);
1128 cl2->SetType(5);
1129 fgPairs[ip2*nn+in] =5;
1130 if ((pos[ip2].GetNd()+neg[in].GetNd())>6){ //multi cluster
1131 cl2->SetType(6);
1132 fgPairs[ip2*nn+in] =6;
1133 }
1134 }
1135 else{
1136 cl2 = new AliITSRecPoint(milab,lp,info);
1137 cl2->SetChargeRatio(ratio);
1138 cl2->SetType(5);
1139 fgPairs[ip2*nn+in] =5;
1140 if ((pos[ip2].GetNd()+neg[in].GetNd())>6){ //multi cluster
1141 cl2->SetType(6);
1142 fgPairs[ip2*nn+in] =6;
1143 }
1144
a86176e3 1145 fDetTypeRec->AddRecPoint(*cl2);
1146 }
1147 ncl++;
a64f9843 1148 }
1149
a86176e3 1150 cused1[ip]++;
1151 cused1[ip2]++;
1152 cused2[in]++;
a64f9843 1153
1154 } // charge matching condition
1155
1156 } // 2 Pside cross 1 Nside
1157 } // loop over Pside clusters
a86176e3 1158
a64f9843 1159
1160
1161 //
1162 for (Int_t jn=0;jn<nn;jn++){
1163 if (cused2[jn]) continue;
8be4e1b1 1164 Float_t xbest=1000,zbest=1000,qbest=0;
a64f9843 1165 // select "silber" cluster
1166 if ( cpositive[jn]==1 && cnegative[positivepair[10*jn]]==2){
1167 Int_t ip = positivepair[10*jn];
1168 Int_t jn2 = negativepair[10*ip];
1169 if (jn2==jn) jn2 = negativepair[10*ip+1];
1170 Float_t pcharge = neg[jn].GetQ()+neg[jn2].GetQ();
04366a57 1171 //
a64f9843 1172
d695268b 1173
1174 ratio = (pcharge-pos[ip].GetQ())/(pcharge+pos[ip].GetQ());
1175 if ( (TMath::Abs(ratio)<0.2) && (pcharge!=0) ) {
1176
1177 /*
2069484c 1178 if ( (TMath::Abs(pcharge-pos[ip].GetQ())<30) && // charge matching
a64f9843 1179 (pcharge!=0) ) { // reject combinations of bad strips
d695268b 1180 */
1181
1182
a64f9843 1183 //
1184 // add first pair
1185 // if (!(cused1[ip]||cused2[jn])){
1186 if ( (fgPairs[ip*nn+jn]==100) && (neg[jn].GetQ()) ) { //
1187
d036ccd3 1188 Float_t yn=neg[jn].GetY();
8be4e1b1 1189 Float_t yp=pos[ip].GetY();
2069484c 1190
8be4e1b1 1191 Float_t xt, zt;
1192 seg->GetPadCxz(yn, yp, xt, zt);
1193
1194 xbest=xt; zbest=zt;
2069484c 1195
a64f9843 1196 qbest =neg[jn].GetQ();
d036ccd3 1197
a64f9843 1198 {
8be4e1b1 1199 Double_t loc[3]={xbest,0.,zbest},trk[3]={0.,0.,0.};
a64f9843 1200 mT2L->MasterToLocal(loc,trk);
1201 lp[0]=trk[1];
1202 lp[1]=trk[2];
b4704be3 1203 }
04366a57 1204
1205 lp[4]=qbest; //Q
1206 for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2;
1207 for (Int_t ilab=0;ilab<3;ilab++){
1208 milab[ilab] = pos[ip].GetLabel(ilab);
1209 milab[ilab+3] = neg[jn].GetLabel(ilab);
1210 }
1211 //
1212 CheckLabels2(milab);
1213 ratio = (pos[ip].GetQ()-neg[jn].GetQ())/(pos[ip].GetQ()+neg[jn].GetQ());
1214 milab[3]=(((ip<<10) + jn)<<10) + idet; // pos|neg|det
1215 Int_t info[3] = {pos[ip].GetNd(),neg[jn].GetNd(),fNlayer[fModule]};
1216
d695268b 1217 lp[2]=0.0022*0.0022; //SigmaY2
1218 lp[3]=0.110*0.110; //SigmaZ2
1219 // out-of-diagonal element of covariance matrix
1220 if( (info[0]==1) && (info[1]==1) ) lp[5]=-0.00012;
1221 else if ( (info[0]>1) && (info[1]>1) ) {
1222 lp[2]=0.0016*0.0016; //SigmaY2
1223 lp[3]=0.08*0.08; //SigmaZ2
1224 lp[5]=-0.00006;
1225 }
1226 else {
7101948c 1227 lp[3]=0.093*0.093;
d695268b 1228 if (info[0]==1) { lp[5]=-0.00014;}
1229 else { lp[2]=0.0017*0.0017; lp[5]=-0.00004;}
1230 }
1231
00a7cc50 1232 AliITSRecPoint * cl2;
1233 if(clusters){
75fb37cc 1234 cl2 = new (cl[ncl]) AliITSRecPoint(milab,lp,info);
308b5ea4 1235
00a7cc50 1236 cl2->SetChargeRatio(ratio);
1237 cl2->SetType(7);
308b5ea4 1238 fgPairs[ip*nn+jn] =7;
00a7cc50 1239 if ((pos[ip].GetNd()+neg[jn].GetNd())>6){ //multi cluster
1240 cl2->SetType(8);
308b5ea4 1241 fgPairs[ip*nn+jn]=8;
00a7cc50 1242 }
1243
1244 }
04366a57 1245 else{
75fb37cc 1246 cl2 = new AliITSRecPoint(milab,lp,info);
00a7cc50 1247 cl2->SetChargeRatio(ratio);
1248 cl2->SetType(7);
308b5ea4 1249 fgPairs[ip*nn+jn] =7;
00a7cc50 1250 if ((pos[ip].GetNd()+neg[jn].GetNd())>6){ //multi cluster
1251 cl2->SetType(8);
308b5ea4 1252 fgPairs[ip*nn+jn]=8;
00a7cc50 1253 }
1254
1255 fDetTypeRec->AddRecPoint(*cl2);
04366a57 1256 }
1257 ncl++;
04366a57 1258 }
1259 //
1260 // add second pair
1261 // if (!(cused1[ip]||cused2[jn2])){
a64f9843 1262 if ( (fgPairs[ip*nn+jn2]==100)&&(neg[jn2].GetQ() ) ) { //
1263
d036ccd3 1264 Float_t yn=neg[jn2].GetY();
1265 Double_t yp=pos[ip].GetY();
2069484c 1266
8be4e1b1 1267 Float_t xt, zt;
1268 seg->GetPadCxz(yn, yp, xt, zt);
1269
1270 xbest=xt; zbest=zt;
2069484c 1271
04366a57 1272 qbest =neg[jn2].GetQ();
d036ccd3 1273
b4704be3 1274 {
8be4e1b1 1275 Double_t loc[3]={xbest,0.,zbest},trk[3]={0.,0.,0.};
b4704be3 1276 mT2L->MasterToLocal(loc,trk);
1277 lp[0]=trk[1];
1278 lp[1]=trk[2];
1279 }
d695268b 1280
04366a57 1281 lp[4]=qbest; //Q
1282 for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2;
1283 for (Int_t ilab=0;ilab<3;ilab++){
1284 milab[ilab] = pos[ip].GetLabel(ilab);
1285 milab[ilab+3] = neg[jn2].GetLabel(ilab);
1286 }
1287 //
1288 CheckLabels2(milab);
1289 ratio = (pos[ip].GetQ()-neg[jn2].GetQ())/(pos[ip].GetQ()+neg[jn2].GetQ());
1290 milab[3]=(((ip<<10) + jn2)<<10) + idet; // pos|neg|det
1291 Int_t info[3] = {pos[ip].GetNd(),neg[jn2].GetNd(),fNlayer[fModule]};
d695268b 1292
1293 lp[2]=0.0022*0.0022; //SigmaY2
1294 lp[3]=0.110*0.110; //SigmaZ2
1295 // out-of-diagonal element of covariance matrix
1296 if( (info[0]==1) && (info[1]==1) ) lp[5]=-0.00012;
1297 else if ( (info[0]>1) && (info[1]>1) ) {
1298 lp[2]=0.0016*0.0016; //SigmaY2
1299 lp[3]=0.08*0.08; //SigmaZ2
1300 lp[5]=-0.00006;
1301 }
1302 else {
7101948c 1303 lp[3]=0.093*0.093;
d695268b 1304 if (info[0]==1) { lp[5]=-0.00014;}
1305 else { lp[2]=0.0017*0.0017; lp[5]=-0.00004;}
1306 }
1307
00a7cc50 1308 AliITSRecPoint * cl2;
1309 if(clusters){
75fb37cc 1310 cl2 = new (cl[ncl]) AliITSRecPoint(milab,lp,info);
308b5ea4 1311
308b5ea4 1312
00a7cc50 1313 cl2->SetChargeRatio(ratio);
308b5ea4 1314 fgPairs[ip*nn+jn2]=7;
00a7cc50 1315 cl2->SetType(7);
1316 if ((pos[ip].GetNd()+neg[jn2].GetNd())>6){ //multi cluster
1317 cl2->SetType(8);
308b5ea4 1318 fgPairs[ip*nn+jn2]=8;
00a7cc50 1319 }
1320
1321 }
04366a57 1322 else{
75fb37cc 1323 cl2 = new AliITSRecPoint(milab,lp,info);
00a7cc50 1324 cl2->SetChargeRatio(ratio);
308b5ea4 1325 fgPairs[ip*nn+jn2]=7;
00a7cc50 1326 cl2->SetType(7);
1327 if ((pos[ip].GetNd()+neg[jn2].GetNd())>6){ //multi cluster
1328 cl2->SetType(8);
308b5ea4 1329 fgPairs[ip*nn+jn2]=8;
00a7cc50 1330 }
308b5ea4 1331
00a7cc50 1332 fDetTypeRec->AddRecPoint(*cl2);
04366a57 1333 }
1334
1335 ncl++;
04366a57 1336 }
1337 cused1[ip]++;
1338 cused2[jn]++;
1339 cused2[jn2]++;
308b5ea4 1340
a64f9843 1341 } // charge matching condition
1342
1343 } // 2 Nside cross 1 Pside
1344 } // loop over Pside clusters
1345
1346
1347
1348 for (Int_t ip=0;ip<np;ip++){
d695268b 1349
1350 if(cused1[ip]) continue;
1351
1352
8be4e1b1 1353 Float_t xbest=1000,zbest=1000,qbest=0;
04366a57 1354 //
a64f9843 1355 // 2x2 clusters
04366a57 1356 //
d695268b 1357 if ( (cnegative[ip]==2) && cpositive[negativepair[10*ip]]==2){
1358 Float_t minchargediff =4.;
1359 Float_t minchargeratio =0.2;
1360
1361 Int_t j=-1;
1362 for (Int_t di=0;di<cnegative[ip];di++){
1363 Int_t jc = negativepair[ip*10+di];
1364 Float_t chargedif = pos[ip].GetQ()-neg[jc].GetQ();
1365 ratio = (pos[ip].GetQ()-neg[jc].GetQ())/(pos[ip].GetQ()+neg[jc].GetQ());
1366 //if (TMath::Abs(chargedif)<minchargediff){
1367 if (TMath::Abs(ratio)<0.2){
1368 j =jc;
1369 minchargediff = TMath::Abs(chargedif);
1370 minchargeratio = TMath::Abs(ratio);
1371 }
1372 }
1373 if (j<0) continue; // not proper cluster
1374
1375
1376 Int_t count =0;
1377 for (Int_t di=0;di<cnegative[ip];di++){
1378 Int_t jc = negativepair[ip*10+di];
1379 Float_t chargedif = pos[ip].GetQ()-neg[jc].GetQ();
1380 if (TMath::Abs(chargedif)<minchargediff+3.) count++;
1381 }
1382 if (count>1) continue; // more than one "proper" cluster for positive
1383 //
1384
1385 count =0;
1386 for (Int_t dj=0;dj<cpositive[j];dj++){
1387 Int_t ic = positivepair[j*10+dj];
1388 Float_t chargedif = pos[ic].GetQ()-neg[j].GetQ();
1389 if (TMath::Abs(chargedif)<minchargediff+3.) count++;
1390 }
1391 if (count>1) continue; // more than one "proper" cluster for negative
1392
1393 Int_t jp = 0;
1394
1395 count =0;
1396 for (Int_t dj=0;dj<cnegative[jp];dj++){
1397 Int_t ic = positivepair[jp*10+dj];
1398 Float_t chargedif = pos[ic].GetQ()-neg[jp].GetQ();
1399 if (TMath::Abs(chargedif)<minchargediff+4.) count++;
1400 }
1401 if (count>1) continue;
1402 if (fgPairs[ip*nn+j]<100) continue;
1403 //
1404
1405
1406
1407 //almost gold clusters
1408 Float_t yp=pos[ip].GetY();
1409 Float_t yn=neg[j].GetY();
1410 Float_t xt, zt;
1411 seg->GetPadCxz(yn, yp, xt, zt);
1412 xbest=xt; zbest=zt;
1413 qbest=0.5*(pos[ip].GetQ()+neg[j].GetQ());
1414 {
1415 Double_t loc[3]={xbest,0.,zbest},trk[3]={0.,0.,0.};
1416 mT2L->MasterToLocal(loc,trk);
1417 lp[0]=trk[1];
1418 lp[1]=trk[2];
1419 }
1420 lp[4]=qbest; //Q
1421 for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2;
1422 for (Int_t ilab=0;ilab<3;ilab++){
1423 milab[ilab] = pos[ip].GetLabel(ilab);
1424 milab[ilab+3] = neg[j].GetLabel(ilab);
1425 }
1426 //
1427 CheckLabels2(milab);
1428 if ((neg[j].GetQ()==0)&&(pos[ip].GetQ()==0)) continue; // reject crosses of bad strips!!
1429 ratio = (pos[ip].GetQ()-neg[j].GetQ())/(pos[ip].GetQ()+neg[j].GetQ());
1430 milab[3]=(((ip<<10) + j)<<10) + idet; // pos|neg|det
1431 Int_t info[3] = {pos[ip].GetNd(),neg[j].GetNd(),fNlayer[fModule]};
1432
1433 lp[2]=0.0022*0.0022; //SigmaY2
1434 lp[3]=0.110*0.110; //SigmaZ2
1435 // out-of-diagonal element of covariance matrix
1436 if( (info[0]==1) && (info[1]==1) ) lp[5]=-0.00012;
1437 else if ( (info[0]>1) && (info[1]>1) ) {
1438 lp[2]=0.0016*0.0016; //SigmaY2
1439 lp[3]=0.08*0.08; //SigmaZ2
1440 lp[5]=-0.00006;
1441 }
1442 else {
7101948c 1443 lp[3]=0.093*0.093;
d695268b 1444 if (info[0]==1) { lp[5]=-0.00014;}
1445 else { lp[2]=0.0017*0.0017; lp[5]=-0.00004;}
1446 }
1447
1448 AliITSRecPoint * cl2;
1449 if(clusters){
1450 cl2 = new (cl[ncl]) AliITSRecPoint(milab,lp,info);
1451
1452 cl2->SetChargeRatio(ratio);
1453 cl2->SetType(10);
1454 fgPairs[ip*nn+j]=10;
1455 if ((pos[ip].GetNd()+neg[j].GetNd())>6){ //multi cluster
1456 cl2->SetType(11);
1457 fgPairs[ip*nn+j]=11;
1458 }
1459 cused1[ip]++;
1460 cused2[j]++;
1461 }
1462 else{
1463 cl2 = new AliITSRecPoint(milab,lp,info);
1464 cl2->SetChargeRatio(ratio);
1465 cl2->SetType(10);
1466 fgPairs[ip*nn+j]=10;
1467 if ((pos[ip].GetNd()+neg[j].GetNd())>6){ //multi cluster
1468 cl2->SetType(11);
1469 fgPairs[ip*nn+j]=11;
1470 }
1471 cused1[ip]++;
1472 cused2[j]++;
1473
1474 fDetTypeRec->AddRecPoint(*cl2);
1475 }
1476 ncl++;
1477
1478 } // 2X2
1479 } // loop over Pside 1Dclusters
1480
1481
1482
1483 for (Int_t ip=0;ip<np;ip++){
1484
1485 if(cused1[ip]) continue;
1486
1487
1488 Float_t xbest=1000,zbest=1000,qbest=0;
1489 //
1490 // manyxmany clusters
1491 //
a64f9843 1492 if ( (cnegative[ip]<5) && cpositive[negativepair[10*ip]]<5){
1493 Float_t minchargediff =4.;
1494 Int_t j=-1;
1495 for (Int_t di=0;di<cnegative[ip];di++){
1496 Int_t jc = negativepair[ip*10+di];
1497 Float_t chargedif = pos[ip].GetQ()-neg[jc].GetQ();
1498 if (TMath::Abs(chargedif)<minchargediff){
1499 j =jc;
1500 minchargediff = TMath::Abs(chargedif);
1501 }
00a7cc50 1502 }
a64f9843 1503 if (j<0) continue; // not proper cluster
1504
1505 Int_t count =0;
1506 for (Int_t di=0;di<cnegative[ip];di++){
1507 Int_t jc = negativepair[ip*10+di];
1508 Float_t chargedif = pos[ip].GetQ()-neg[jc].GetQ();
1509 if (TMath::Abs(chargedif)<minchargediff+3.) count++;
00a7cc50 1510 }
a64f9843 1511 if (count>1) continue; // more than one "proper" cluster for positive
1512 //
00a7cc50 1513
a64f9843 1514 count =0;
1515 for (Int_t dj=0;dj<cpositive[j];dj++){
1516 Int_t ic = positivepair[j*10+dj];
1517 Float_t chargedif = pos[ic].GetQ()-neg[j].GetQ();
1518 if (TMath::Abs(chargedif)<minchargediff+3.) count++;
1519 }
1520 if (count>1) continue; // more than one "proper" cluster for negative
1521
1522 Int_t jp = 0;
1523
1524 count =0;
1525 for (Int_t dj=0;dj<cnegative[jp];dj++){
1526 Int_t ic = positivepair[jp*10+dj];
1527 Float_t chargedif = pos[ic].GetQ()-neg[jp].GetQ();
1528 if (TMath::Abs(chargedif)<minchargediff+4.) count++;
1529 }
1530 if (count>1) continue;
1531 if (fgPairs[ip*nn+j]<100) continue;
1532 //
1533
1534 //almost gold clusters
d036ccd3 1535 Float_t yp=pos[ip].GetY();
1536 Float_t yn=neg[j].GetY();
1537
8be4e1b1 1538
1539 Float_t xt, zt;
1540 seg->GetPadCxz(yn, yp, xt, zt);
1541
1542 xbest=xt; zbest=zt;
1543
a64f9843 1544 qbest=0.5*(pos[ip].GetQ()+neg[j].GetQ());
d036ccd3 1545
a64f9843 1546 {
8be4e1b1 1547 Double_t loc[3]={xbest,0.,zbest},trk[3]={0.,0.,0.};
a64f9843 1548 mT2L->MasterToLocal(loc,trk);
1549 lp[0]=trk[1];
1550 lp[1]=trk[2];
1551 }
a64f9843 1552 lp[4]=qbest; //Q
1553 for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2;
1554 for (Int_t ilab=0;ilab<3;ilab++){
1555 milab[ilab] = pos[ip].GetLabel(ilab);
1556 milab[ilab+3] = neg[j].GetLabel(ilab);
1557 }
1558 //
1559 CheckLabels2(milab);
db6e54cd 1560 if ((neg[j].GetQ()==0)&&(pos[ip].GetQ()==0)) continue; // reject crosses of bad strips!!
a64f9843 1561 ratio = (pos[ip].GetQ()-neg[j].GetQ())/(pos[ip].GetQ()+neg[j].GetQ());
1562 milab[3]=(((ip<<10) + j)<<10) + idet; // pos|neg|det
1563 Int_t info[3] = {pos[ip].GetNd(),neg[j].GetNd(),fNlayer[fModule]};
d695268b 1564
1565 lp[2]=0.0022*0.0022; //SigmaY2
1566 lp[3]=0.110*0.110; //SigmaZ2
1567 // out-of-diagonal element of covariance matrix
1568 if( (info[0]==1) && (info[1]==1) ) lp[5]=-0.00012;
1569 else if ( (info[0]>1) && (info[1]>1) ) {
1570 lp[2]=0.0016*0.0016; //SigmaY2
1571 lp[3]=0.08*0.08; //SigmaZ2
1572 lp[5]=-0.00006;
1573 }
1574 else {
7101948c 1575 lp[3]=0.093*0.093;
d695268b 1576 if (info[0]==1) { lp[5]=-0.00014;}
1577 else { lp[2]=0.0017*0.0017; lp[5]=-0.00004;}
1578 }
1579
a64f9843 1580 AliITSRecPoint * cl2;
1581 if(clusters){
1582 cl2 = new (cl[ncl]) AliITSRecPoint(milab,lp,info);
d036ccd3 1583
a64f9843 1584 cl2->SetChargeRatio(ratio);
d695268b 1585 cl2->SetType(12);
1586 fgPairs[ip*nn+j]=12;
a64f9843 1587 if ((pos[ip].GetNd()+neg[j].GetNd())>6){ //multi cluster
d695268b 1588 cl2->SetType(13);
1589 fgPairs[ip*nn+j]=13;
a64f9843 1590 }
1591 cused1[ip]++;
1592 cused2[j]++;
1593 }
1594 else{
1595 cl2 = new AliITSRecPoint(milab,lp,info);
1596 cl2->SetChargeRatio(ratio);
d695268b 1597 cl2->SetType(12);
1598 fgPairs[ip*nn+j]=12;
a64f9843 1599 if ((pos[ip].GetNd()+neg[j].GetNd())>6){ //multi cluster
d695268b 1600 cl2->SetType(13);
1601 fgPairs[ip*nn+j]=13;
a64f9843 1602 }
1603 cused1[ip]++;
1604 cused2[j]++;
1605
a64f9843 1606 fDetTypeRec->AddRecPoint(*cl2);
1607 }
1608 ncl++;
1609
1610 } // manyXmany
1611 } // loop over Pside 1Dclusters
1612
a64f9843 1613 } // use charge matching
1614
04366a57 1615
a64f9843 1616 // recover all the other crosses
04366a57 1617 //
1618 for (Int_t i=0; i<np; i++) {
8be4e1b1 1619 Float_t xbest=1000,zbest=1000,qbest=0;
d036ccd3 1620 Float_t yp=pos[i].GetY();
a64f9843 1621 if ((pos[i].GetQ()>0)&&(pos[i].GetQ()<3)) continue;
04366a57 1622 for (Int_t j=0; j<nn; j++) {
1623 // for (Int_t di = 0;di<cpositive[i];di++){
1624 // Int_t j = negativepair[10*i+di];
a64f9843 1625 if ((neg[j].GetQ()>0)&&(neg[j].GetQ()<3)) continue;
1626
1627 if ((neg[j].GetQ()==0)&&(pos[i].GetQ()==0)) continue; // reject crosses of bad strips!!
1628
04366a57 1629 if (cused2[j]||cused1[i]) continue;
308b5ea4 1630 if (fgPairs[i*nn+j]>0 &&fgPairs[i*nn+j]<100) continue;
04366a57 1631 ratio = (pos[i].GetQ()-neg[j].GetQ())/(pos[i].GetQ()+neg[j].GetQ());
d036ccd3 1632 Float_t yn=neg[j].GetY();
1633
8be4e1b1 1634 Float_t xt, zt;
1635 seg->GetPadCxz(yn, yp, xt, zt);
d036ccd3 1636
8be4e1b1 1637 if (TMath::Abs(xt)<hwSSD+0.01)
1638 if (TMath::Abs(zt)<hlSSD+0.01*(neg[j].GetNd()+pos[i].GetNd())) {
1639 xbest=xt; zbest=zt;
1640
04366a57 1641 qbest=0.5*(pos[i].GetQ()+neg[j].GetQ());
d036ccd3 1642
b4704be3 1643 {
8be4e1b1 1644 Double_t loc[3]={xbest,0.,zbest},trk[3]={0.,0.,0.};
b4704be3 1645 mT2L->MasterToLocal(loc,trk);
1646 lp[0]=trk[1];
1647 lp[1]=trk[2];
1648 }
04366a57 1649 lp[4]=qbest; //Q
1650 for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2;
1651 for (Int_t ilab=0;ilab<3;ilab++){
1652 milab[ilab] = pos[i].GetLabel(ilab);
1653 milab[ilab+3] = neg[j].GetLabel(ilab);
1654 }
1655 //
1656 CheckLabels2(milab);
1657 milab[3]=(((i<<10) + j)<<10) + idet; // pos|neg|det
1658 Int_t info[3] = {pos[i].GetNd(),neg[j].GetNd(),fNlayer[fModule]};
d695268b 1659
1660 lp[2]=0.0022*0.0022; //SigmaY2
1661 lp[3]=0.110*0.110; //SigmaZ2
1662 // out-of-diagonal element of covariance matrix
1663 if( (info[0]==1) && (info[1]==1) ) lp[5]=-0.00012;
1664 else if ( (info[0]>1) && (info[1]>1) ) {
1665 lp[2]=0.0016*0.0016; //SigmaY2
1666 lp[3]=0.08*0.08; //SigmaZ2
1667 lp[5]=-0.00006;
1668 }
1669 else {
7101948c 1670 lp[3]=0.093*0.093;
d695268b 1671 if (info[0]==1) { lp[5]=-0.00014;}
1672 else { lp[2]=0.0017*0.0017; lp[5]=-0.00004;}
1673 }
1674
00a7cc50 1675 AliITSRecPoint * cl2;
1676 if(clusters){
75fb37cc 1677 cl2 = new (cl[ncl]) AliITSRecPoint(milab,lp,info);
308b5ea4 1678
00a7cc50 1679 cl2->SetChargeRatio(ratio);
1680 cl2->SetType(100+cpositive[j]+cnegative[i]);
a64f9843 1681
1682 if(pos[i].GetQ()==0) cl2->SetType(200+cpositive[j]+cnegative[i]);
1683 if(neg[j].GetQ()==0) cl2->SetType(300+cpositive[j]+cnegative[i]);
1684
00a7cc50 1685 }
04366a57 1686 else{
75fb37cc 1687 cl2 = new AliITSRecPoint(milab,lp,info);
00a7cc50 1688 cl2->SetChargeRatio(ratio);
1689 cl2->SetType(100+cpositive[j]+cnegative[i]);
308b5ea4 1690
a64f9843 1691 if(pos[i].GetQ()==0) cl2->SetType(200+cpositive[j]+cnegative[i]);
1692 if(neg[j].GetQ()==0) cl2->SetType(300+cpositive[j]+cnegative[i]);
308b5ea4 1693
00a7cc50 1694 fDetTypeRec->AddRecPoint(*cl2);
04366a57 1695 }
1696 ncl++;
04366a57 1697 }
1698 }
1699 }
d695268b 1700
1701
d695268b 1702
7101948c 1703 if(repa->GetUseBadChannelsInClusterFinderSSD()==kTRUE) {
1704
1705 //---------------------------------------------------------
1706 // recover crosses of good 1D clusters with bad strips on the other side
1707 // Note1: at first iteration skip modules with a bad side (or almost), (would produce too many fake!)
1708 // Note2: for modules with a bad side see below
1709
1710 AliITSCalibrationSSD* cal = (AliITSCalibrationSSD*)GetResp(fModule);
1711 Int_t countPbad=0, countNbad=0;
1712 for(Int_t ib=0; ib<768; ib++) {
1713 if(cal->IsPChannelBad(ib)) countPbad++;
1714 if(cal->IsNChannelBad(ib)) countNbad++;
1715 }
1716 // AliInfo(Form("module %d has %d P- and %d N-bad strips",fModule,countPbad,countNbad));
d695268b 1717
7101948c 1718 if( (countPbad<100) && (countNbad<100) ) { // no bad side!!
d695268b 1719
7101948c 1720 for (Int_t i=0; i<np; i++) { // loop over Nside 1Dclusters with no crosses
1721 if(cnegative[i]) continue; // if intersecting Pside clusters continue;
d695268b 1722
7101948c 1723 // for(Int_t ib=0; ib<768; ib++) { // loop over all Pstrips
1724 for(Int_t ib=15; ib<753; ib++) { // loop over all Pstrips
d695268b 1725
7101948c 1726 if(cal->IsPChannelBad(ib)) { // check if strips is bad
1727 Float_t yN=pos[i].GetY();
1728 Float_t xt, zt;
1729 seg->GetPadCxz(1.*ib, yN, xt, zt);
d695268b 1730
7101948c 1731 //----------
1732 // bad Pstrip is crossing the Nside 1Dcluster -> create recpoint
1733 //
1734 if ( (TMath::Abs(xt)<hwSSD+0.01) && (TMath::Abs(zt)<hlSSD+0.01) ) {
1735 Double_t loc[3]={xt,0.,zt},trk[3]={0.,0.,0.};
1736 mT2L->MasterToLocal(loc,trk);
1737 lp[0]=trk[1];
1738 lp[1]=trk[2];
1739 lp[4]=pos[i].GetQ(); //Q
1740 for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2;
1741 for (Int_t ilab=0;ilab<3;ilab++) milab[ilab] = pos[i].GetLabel(ilab);
1742 CheckLabels2(milab);
1743 milab[3]=( (i<<10) << 10 ) + idet; // pos|neg|det
1744 Int_t info[3] = {pos[i].GetNd(),0,fNlayer[fModule]};
1745
1746 // out-of-diagonal element of covariance matrix
1747 if (info[0]==1) lp[5]=0.0065;
1748 else lp[5]=0.0093;
1749
1750 lp[2]=0.0022*0.0022; //SigmaY2
1751 lp[3]=0.110*0.110; //SigmaZ2
1752 lp[5]=-0.00012; // out-of-diagonal element of covariance matrix
1753
1754 AliITSRecPoint * cl2;
1755 if(clusters){
1756 cl2 = new (cl[ncl]) AliITSRecPoint(milab,lp,info);
1757 cl2->SetChargeRatio(1.);
d695268b 1758 cl2->SetType(50);
7101948c 1759 }
1760 else{
1761 cl2 = new AliITSRecPoint(milab,lp,info);
1762 cl2->SetChargeRatio(1.);
1763 cl2->SetType(50);
1764 fDetTypeRec->AddRecPoint(*cl2);
1765 }
1766 ncl++;
1767 } // cross is within the detector
1768 //
1769 //--------------
1770
1771 } // bad Pstrip
d695268b 1772
7101948c 1773 } // end loop over Pstrips
d695268b 1774
7101948c 1775 } // end loop over Nside 1D clusters
d695268b 1776
7101948c 1777 for (Int_t j=0; j<nn; j++) { // loop over Pside 1D clusters with no crosses
1778 if(cpositive[j]) continue;
1779
1780 // for(Int_t ib=0; ib<768; ib++) { // loop over all Nside strips
1781 for(Int_t ib=15; ib<753; ib++) { // loop over all Nside strips
1782
1783 if(cal->IsNChannelBad(ib)) { // check if strip is bad
1784 Float_t yP=neg[j].GetY();
1785 Float_t xt, zt;
1786 seg->GetPadCxz(yP, 1.*ib, xt, zt);
1787
1788 //----------
1789 // bad Nstrip is crossing the Pside 1Dcluster -> create recpoint
1790 //
1791 if ( (TMath::Abs(xt)<hwSSD+0.01) && (TMath::Abs(zt)<hlSSD+0.01) ) {
1792 Double_t loc[3]={xt,0.,zt},trk[3]={0.,0.,0.};
1793 mT2L->MasterToLocal(loc,trk);
1794 lp[0]=trk[1];
1795 lp[1]=trk[2];
1796 lp[4]=neg[j].GetQ(); //Q
1797 for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2;
1798 for (Int_t ilab=0;ilab<3;ilab++) milab[ilab] = neg[j].GetLabel(ilab);
1799 CheckLabels2(milab);
1800 milab[3]=( j << 10 ) + idet; // pos|neg|det
1801 Int_t info[3]={0,(Int_t)neg[j].GetNd(),fNlayer[fModule]};
1802
1803 lp[2]=0.0022*0.0022; //SigmaY2
1804 lp[3]=0.110*0.110; //SigmaZ2
1805 lp[5]=-0.00012; // out-of-diagonal element of covariance matrix
1806
1807 AliITSRecPoint * cl2;
1808 if(clusters){
1809 cl2 = new (cl[ncl]) AliITSRecPoint(milab,lp,info);
1810 cl2->SetChargeRatio(1.);
1811 cl2->SetType(60);
1812 }
1813 else{
1814 cl2 = new AliITSRecPoint(milab,lp,info);
1815 cl2->SetChargeRatio(1.);
1816 cl2->SetType(60);
1817 fDetTypeRec->AddRecPoint(*cl2);
1818 }
1819 ncl++;
1820 } // cross is within the detector
1821 //
1822 //--------------
1823
1824 } // bad Nstrip
1825 } // end loop over Nstrips
1826 } // end loop over Pside 1D clusters
1827
1828 } // no bad sides
1829
1830 //---------------------------------------------------------
d695268b 1831
7101948c 1832 else if( (countPbad>700) && (countNbad<100) ) { // bad Pside!!
d695268b 1833
7101948c 1834 for (Int_t i=0; i<np; i++) { // loop over Nside 1Dclusters with no crosses
1835 if(cnegative[i]) continue; // if intersecting Pside clusters continue;
d695268b 1836
7101948c 1837 Float_t xt, zt;
1838 Float_t yN=pos[i].GetY();
1839 Float_t yP=0.;
1840 if (seg->GetLayer()==5) yP = yN + (7.6/1.9);
1841 else yP = yN - (7.6/1.9);
1842 seg->GetPadCxz(yP, yN, xt, zt);
1843
1844 if ( (TMath::Abs(xt)<hwSSD+0.01) && (TMath::Abs(zt)<hlSSD+0.01) ) {
1845 Double_t loc[3]={xt,0.,zt},trk[3]={0.,0.,0.};
1846 mT2L->MasterToLocal(loc,trk);
1847 lp[0]=trk[1];
1848 lp[1]=trk[2];
1849 lp[4]=pos[i].GetQ(); //Q
1850 for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2;
1851 for (Int_t ilab=0;ilab<3;ilab++) milab[ilab] = pos[i].GetLabel(ilab);
1852 CheckLabels2(milab);
1853 milab[3]=( (i<<10) << 10 ) + idet; // pos|neg|det
1854 Int_t info[3] = {(Int_t)pos[i].GetNd(),0,fNlayer[fModule]};
d695268b 1855
7101948c 1856 lp[2]=0.031*0.031; //SigmaY2
1857 lp[3]=1.15*1.15; //SigmaZ2
1858 lp[5]=-0.036;
d695268b 1859
7101948c 1860 AliITSRecPoint * cl2;
1861 if(clusters){
1862 cl2 = new (cl[ncl]) AliITSRecPoint(milab,lp,info);
1863 cl2->SetChargeRatio(1.);
1864 cl2->SetType(70);
1865 }
1866 else{
1867 cl2 = new AliITSRecPoint(milab,lp,info);
1868 cl2->SetChargeRatio(1.);
1869 cl2->SetType(70);
1870 fDetTypeRec->AddRecPoint(*cl2);
1871 }
1872 ncl++;
1873 } // cross is within the detector
1874 //
1875 //--------------
1876
1877 } // end loop over Nside 1D clusters
1878
1879 } // bad Pside module
d695268b 1880
7101948c 1881 else if( (countNbad>700) && (countPbad<100) ) { // bad Nside!!
d695268b 1882
7101948c 1883 for (Int_t j=0; j<nn; j++) { // loop over Pside 1D clusters with no crosses
1884 if(cpositive[j]) continue;
1885
1886 Float_t xt, zt;
1887 Float_t yP=neg[j].GetY();
1888 Float_t yN=0.;
1889 if (seg->GetLayer()==5) yN = yP - (7.6/1.9);
1890 else yN = yP + (7.6/1.9);
1891 seg->GetPadCxz(yP, yN, xt, zt);
1892
1893 if ( (TMath::Abs(xt)<hwSSD+0.01) && (TMath::Abs(zt)<hlSSD+0.01) ) {
1894 Double_t loc[3]={xt,0.,zt},trk[3]={0.,0.,0.};
1895 mT2L->MasterToLocal(loc,trk);
1896 lp[0]=trk[1];
1897 lp[1]=trk[2];
1898 lp[4]=neg[j].GetQ(); //Q
1899 for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2;
1900 for (Int_t ilab=0;ilab<3;ilab++) milab[ilab] = neg[j].GetLabel(ilab);
1901 CheckLabels2(milab);
1902 milab[3]=( j << 10 ) + idet; // pos|neg|det
1903 Int_t info[3] = {0,(Int_t)neg[j].GetNd(),fNlayer[fModule]};
d695268b 1904
7101948c 1905 lp[2]=0.0085*0.0085; //SigmaY2
1906 lp[3]=1.15*1.15; //SigmaZ2
1907 lp[5]=0.0093;
1908
1909 AliITSRecPoint * cl2;
1910 if(clusters){
1911 cl2 = new (cl[ncl]) AliITSRecPoint(milab,lp,info);
1912 cl2->SetChargeRatio(1.);
1913 cl2->SetType(80);
1914 }
1915 else{
1916 cl2 = new AliITSRecPoint(milab,lp,info);
1917 cl2->SetChargeRatio(1.);
1918 cl2->SetType(80);
1919 fDetTypeRec->AddRecPoint(*cl2);
1920 }
1921 ncl++;
1922 } // cross is within the detector
1923 //
1924 //--------------
1925
1926 } // end loop over Pside 1D clusters
d695268b 1927
7101948c 1928 } // bad Nside module
d695268b 1929
7101948c 1930 //---------------------------------------------------------
d695268b 1931
7101948c 1932 } // use bad channels
d695268b 1933
d695268b 1934 //cout<<ncl<<" clusters for this module"<<endl;
1935
8be4e1b1 1936 delete [] cnegative;
1937 delete [] cused1;
1938 delete [] negativepair;
1939 delete [] cpositive;
1940 delete [] cused2;
1941 delete [] positivepair;
04366a57 1942
0a56760a 1943}