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