]> git.uio.no Git - u/mrichter/AliRoot.git/blame - ITS/ITSbase/AliITSClusterFinderV2SSD.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Overlaps corrected, new shape of sectors
[u/mrichter/AliRoot.git] / ITS / ITSbase / 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;
15be6f41 188
66b89079 189 Int_t layer = 4;
190 if (fModule>fLastSSD1)
191 layer = 5;
15be6f41 192 */
d695268b 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
15be6f41 513 /* Int_t layer = 4;
66b89079 514 if (fModule>fLastSSD1)
515 layer = 5;
15be6f41 516 */
a64f9843 517
138df073 518 AliITSCalibrationSSD* cal = (AliITSCalibrationSSD*)fDetTypeRec->GetCalibrationModel(fModule);
519 if( !cal ){
520 AliWarning(Form("ITSClustersFinderSSD: No calibration found for module (ddl %d ad %d adc %d)",ddl,ad,adc));
521 continue;
522 }
308b5ea4 523
138df073 524 Float_t dStrip = 0;
525
526 if( repa->GetUseCosmicRunShiftsSSD()) { // Special condition for 2007/2008 cosmic data
527 dStrip = fgkCosmic2008StripShifts[ddl][ad-1];
528 if (TMath::Abs(dStrip) > 1.5){
529 AliWarning(Form("Indexing error in Cosmic calibration: ddl = %d, dStrip %f\n",ddl,dStrip));
530 dStrip = 0;
531 }
308b5ea4 532 }
b42cfa25 533
138df073 534 for( int side=0; side<=1; side++ ){
535
536 Int_t lab[3]={-2,-2,-2};
537 Float_t q = 0.;
538 Float_t y = 0.;
539 Int_t nDigits = 0;
540 Int_t ostrip = -2;
541 Bool_t snFlag = 0;
66b89079 542
543 Float_t dLorentz = 0;
544 if (side==0) { // P-side is neg clust
545 dLorentz = fLorentzShiftN;
546 }
547 else { // N-side is pos clust
548 dLorentz = fLorentzShiftP;
549 }
b42cfa25 550
138df073 551 Int_t n = nStrips[adc][side];
552 for( int istr = 0; istr<n+1; istr++ ){
b42cfa25 553
138df073 554 bool stripOK = 1;
6c8e94cf 555 Int_t strip=0, rwID = 0;
138df073 556 Float_t signal=0.0, noise=0.0, gain=0.0;
308b5ea4 557
138df073 558 if( istr<n ){
559 strip = strips[adc][side][istr][0];
560 signal = strips[adc][side][istr][1];
6c8e94cf 561 rwID = strips[adc][side][istr][2]; // RS
138df073 562 //cout<<"strip "<<adc<<" / "<<side<<": "<<strip<<endl;
563
564 if( cal ){
565 noise = side ?cal->GetNoiseN(strip) :cal->GetNoiseP(strip);
566 gain = side ?cal->GetGainN(strip) :cal->GetGainP(strip);
567 stripOK = ( noise>=1. && signal>=3.0*noise
568 //&& !cal->IsPChannelBad(strip)
569 );
570 }
571 } else stripOK = 0; // end of data
572
573 bool newCluster = ( (abs(strip-ostrip)!=1) || !stripOK );
574
575 if( newCluster ){
576
577 //* Store the previous cluster
578
579 if( nDigits>0 && q>0 && snFlag ){
580
581 if (nClusters1D[side] >= kMaxADCClusters-1 ) {
582 AliWarning("HLT ClustersFinderSSD: Too many 1D clusters !");
583 }else {
584
585 Ali1Dcluster &cluster = clusters1D[side][nClusters1D[side]++];
66b89079 586 cluster.SetY( y / q + dStrip + dLorentz);
138df073 587 cluster.SetQ(q);
588 cluster.SetNd(nDigits);
589 cluster.SetLabels(lab);
590 //cout<<"cluster 1D side "<<side<<": y= "<<y<<" q= "<<q<<" d="<<dStrip<<" Y="<<cluster.GetY()<<endl;
591 //Split suspiciously big cluster
592
593 if( repa->GetUseUnfoldingInClusterFinderSSD()
594 && nDigits > 4 && nDigits < 25
595 ){
66b89079 596 cluster.SetY(y/q + dStrip - 0.25*nDigits + dLorentz);
138df073 597 cluster.SetQ(0.5*q);
598 Ali1Dcluster& cluster2 = clusters1D[side][nClusters1D[side]++];
66b89079 599 cluster2.SetY(y/q + dStrip + 0.25*nDigits + dLorentz);
138df073 600 cluster2.SetQ(0.5*q);
601 cluster2.SetNd(nDigits);
602 cluster2.SetLabels(lab);
603 } // unfolding is on
a86176e3 604 }
308b5ea4 605 }
138df073 606 y = q = 0.;
607 nDigits = 0;
608 snFlag = 0;
609
610 } //* End store the previous cluster
611
612 if( stripOK ){ // add new signal to the cluster
613// signal = (Int_t) (signal * gain); // signal is corrected for gain
614 if( signal>fgkThreshold*noise) snFlag = 1;
615 signal = signal * gain; // signal is corrected for gain
616// if( cal ) signal = (Int_t) cal->ADCToKeV( signal ); // signal is converted in KeV
617 if( cal ) signal = cal->ADCToKeV( signal ); // signal is converted in KeV
618 q += signal; // add digit to current cluster
619 y += strip * signal;
620 nDigits++;
621 //nstat[side]++;
622 ostrip = strip;
6c8e94cf 623 if (fRawID2ClusID) fRawIDRef[side].AddReference(nClusters1D[side],rwID);
a64f9843 624
308b5ea4 625 }
138df073 626 } //* end loop over strips
b42cfa25 627
138df073 628 } //* end loop over ADC sides
308b5ea4 629
a64f9843 630
138df073 631 //* 2D clusterfinder
632 if( nClusters1D[0] && nClusters1D[1] && fModule>=0 ){
633 TClonesArray* clusters = rpc->UncheckedGetClusters(fModule);
634 FindClustersSSD( clusters1D[0], nClusters1D[0], clusters1D[1], nClusters1D[1], clusters);
635 Int_t nClustersn = clusters->GetEntriesFast();
636 nClustersSSD += nClustersn;
308b5ea4 637 }
638
138df073 639 //cout<<"SG: "<<ddl<<" "<<ad<<" "<<adc<<": strips "<<nstat[0]<<"+"<<nstat[1]<<", clusters 1D= "<<nClusters1D[0]<<" + "<<nClusters1D[1]<<", 2D= "<<clusters.size()<<endl;
a64f9843 640
138df073 641 }//* end loop over adc
642
643 }//* end of reconstruction of previous block of 12 modules
644
645 if( newModule ){
646
647 //*
648 //* Clean up arrays and set new module
649 //*
650
651 for( int i=0; i<kNADC; i++ ){
652 nStrips[i][0] = 0;
653 nStrips[i][1] = 0;
654 }
655 ddl = newDDL;
656 ad = newAD;
657 }
658
a64f9843 659
138df073 660 //*
661 //* Exit main loop when there is no more input
662 //*
04366a57 663
138df073 664 if( !next ) break;
665
666 //*
667 //* Fill the current strip information
668 //*
669
670 Int_t adc = input->GetADC();
671 if( adc<0 || adc>=kNADC+2 || (adc>5&&adc<8) ){
672 AliWarning(Form("HLT ClustersFinderSSD: Corrupted data: wrong adc number (%d)", adc));
673 continue;
674 }
a64f9843 675
138df073 676 if( adc>7 ) adc-= 2; // shift ADC numbers 8-13 to 6-11
677
678 Bool_t side = input->GetSideFlag();
679 Int_t strip = input->GetStrip();
680 Int_t signal = input->GetSignal();
681
a86176e3 682
138df073 683 //cout<<"SSD: "<<ddl<<" "<<ad<<" "<<adc<<" "<<side<<" "<<strip<<" : "<<signal<<endl;
04366a57 684
138df073 685 if( strip>767 ){
686 AliWarning(Form("HLT ClustersFinderSSD: Corrupted data: wrong strip number (ddl %d ad %d adc %d side %d, strip %d",
687 ddl, ad, adc, side,strip) );
688 continue;
689 }
690 if (strip < 0) continue;
691
692 int &n = nStrips[adc][side];
693 if( n >0 ){
694 Int_t oldStrip = strips[adc][side][n-1][0];
695
696 if( strip==oldStrip ){
697 AliWarning(Form("HLT ClustersFinderSSD: Corrupted data: duplicated signal: ddl %d ad %d adc %d, side %d, strip %d",
698 ddl, ad, adc, side, strip ));
699 continue;
700 }
701 }
702 strips[adc][side][n][0] = strip;
703 strips[adc][side][n][1] = signal;
6c8e94cf 704 strips[adc][side][n][2] = countRW;
138df073 705 n++;
3a4139a2 706
138df073 707 //cout<<"SSD: "<<input->GetDDL()<<" "<<input->GetAD()<<" "
708 //<<input->GetADC()<<" "<<input->GetSideFlag()<<" "<<((int)input->GetStrip())<<" "<<strip<<" : "<<input->GetSignal()<<endl;
6c8e94cf 709 //
710 countRW++; //RS
138df073 711 } //* End main loop over the input
308b5ea4 712
3b1d8321 713 AliDebug(1,Form("found clusters in ITS SSD: %d", nClustersSSD));
04366a57 714}
715
138df073 716
04366a57 717void AliITSClusterFinderV2SSD::
18f63405 718FindClustersSSD(const Ali1Dcluster* neg, Int_t nn,
719 const Ali1Dcluster* pos, Int_t np,
d1781ed9 720 TClonesArray *clusters) {
04366a57 721 //------------------------------------------------------------
722 // Actual SSD cluster finder
723 //------------------------------------------------------------
b4704be3 724
725 const TGeoHMatrix *mT2L=AliITSgeomTGeo::GetTracking2LocalMatrix(fModule);
726
7101948c 727 //---------------------------------------
728 // load recoparam
729 //
730 static AliITSRecoParam *repa = NULL;
731 if(!repa){
732 repa = (AliITSRecoParam*) AliITSReconstructor::GetRecoParam();
733 if(!repa){
734 repa = AliITSRecoParam::GetHighFluxParam();
735 AliWarning("Using default AliITSRecoParam class");
736 }
737 }
738
138df073 739// TClonesArray &cl=*clusters;
d036ccd3 740
443e42aa 741 AliITSsegmentationSSD *seg = static_cast<AliITSsegmentationSSD*>(fDetTypeRec->GetSegmentationModel(2));
8be4e1b1 742 if (fModule>fLastSSD1)
743 seg->SetLayer(6);
744 else
745 seg->SetLayer(5);
746
747 Float_t hwSSD = seg->Dx()*1e-4/2;
748 Float_t hlSSD = seg->Dz()*1e-4/2;
749
04366a57 750 Int_t idet=fNdet[fModule];
751 Int_t ncl=0;
d036ccd3 752
04366a57 753 //
8be4e1b1 754 Int_t *cnegative = new Int_t[np];
755 Int_t *cused1 = new Int_t[np];
756 Int_t *negativepair = new Int_t[10*np];
757 Int_t *cpositive = new Int_t[nn];
758 Int_t *cused2 = new Int_t[nn];
759 Int_t *positivepair = new Int_t[10*nn];
760 for (Int_t i=0;i<np;i++) {cnegative[i]=0; cused1[i]=0;}
761 for (Int_t i=0;i<nn;i++) {cpositive[i]=0; cused2[i]=0;}
762 for (Int_t i=0;i<10*np;i++) {negativepair[i]=0;}
763 for (Int_t i=0;i<10*nn;i++) {positivepair[i]=0;}
c157c94e 764
308b5ea4 765 if ((np*nn) > fgPairsSize) {
d036ccd3 766
d1781ed9 767 delete [] fgPairs;
d9f9d128 768 fgPairsSize = 2*np*nn;
308b5ea4 769 fgPairs = new Short_t[fgPairsSize];
c157c94e 770 }
308b5ea4 771 memset(fgPairs,0,sizeof(Short_t)*np*nn);
772
04366a57 773 //
774 // find available pairs
775 //
138df073 776 Int_t ncross = 0;
04366a57 777 for (Int_t i=0; i<np; i++) {
d036ccd3 778 Float_t yp=pos[i].GetY();
a64f9843 779 if ( (pos[i].GetQ()>0) && (pos[i].GetQ()<3) ) continue;
04366a57 780 for (Int_t j=0; j<nn; j++) {
a64f9843 781 if ( (neg[j].GetQ()>0) && (neg[j].GetQ()<3) ) continue;
d036ccd3 782 Float_t yn=neg[j].GetY();
783
8be4e1b1 784 Float_t xt, zt;
785 seg->GetPadCxz(yn, yp, xt, zt);
d695268b 786 //cout<<yn<<" "<<yp<<" "<<xt<<" "<<zt<<endl;
d036ccd3 787
5846520e 788 if (TMath::Abs(xt)<hwSSD)
789 if (TMath::Abs(zt)<hlSSD) {
710f576f 790 Int_t in = i*10+cnegative[i];
791 Int_t ip = j*10+cpositive[j];
792 if ((in < 10*np) && (ip < 10*nn)) {
793 negativepair[in] =j; //index
794 positivepair[ip] =i;
795 cnegative[i]++; //counters
138df073 796 cpositive[j]++;
797 ncross++;
710f576f 798 fgPairs[i*nn+j]=100;
799 }
800 else
801 AliError(Form("Index out of range: ip=%d, in=%d",ip,in));
04366a57 802 }
803 }
804 }
308b5ea4 805
eb595a70 806 if (!ncross) {
807 delete [] cnegative;
808 delete [] cused1;
809 delete [] negativepair;
810 delete [] cpositive;
811 delete [] cused2;
812 delete [] positivepair;
813 return;
814 }
138df073 815//why not to allocate memorey here? if(!clusters) clusters = new TClonesArray("AliITSRecPoint", ncross);
816
5846520e 817 /* //
308b5ea4 818 // try to recover points out of but close to the module boundaries
04366a57 819 //
820 for (Int_t i=0; i<np; i++) {
d036ccd3 821 Float_t yp=pos[i].GetY();
a64f9843 822 if ( (pos[i].GetQ()>0) && (pos[i].GetQ()<3) ) continue;
04366a57 823 for (Int_t j=0; j<nn; j++) {
a64f9843 824 if ( (neg[j].GetQ()>0) && (neg[j].GetQ()<3) ) continue;
308b5ea4 825 // if both 1Dclusters have an other cross continue
04366a57 826 if (cpositive[j]&&cnegative[i]) continue;
8be4e1b1 827 Float_t yn=neg[j].GetY();
d036ccd3 828
8be4e1b1 829 Float_t xt, zt;
830 seg->GetPadCxz(yn, yp, xt, zt);
d036ccd3 831
8be4e1b1 832 if (TMath::Abs(xt)<hwSSD+0.1)
833 if (TMath::Abs(zt)<hlSSD+0.15) {
308b5ea4 834 // tag 1Dcluster (eventually will produce low quality recpoint)
04366a57 835 if (cnegative[i]==0) pos[i].SetNd(100); // not available pair
836 if (cpositive[j]==0) neg[j].SetNd(100); // not available pair
710f576f 837 Int_t in = i*10+cnegative[i];
838 Int_t ip = j*10+cpositive[j];
839 if ((in < 10*np) && (ip < 10*nn)) {
840 negativepair[in] =j; //index
841 positivepair[ip] =i;
842 cnegative[i]++; //counters
843 cpositive[j]++;
844 fgPairs[i*nn+j]=100;
845 }
846 else
847 AliError(Form("Index out of range: ip=%d, in=%d",ip,in));
04366a57 848 }
849 }
850 }
5846520e 851 */
308b5ea4 852
04366a57 853 //
d695268b 854 Float_t lp[6];
04366a57 855 Int_t milab[10];
856 Double_t ratio;
857
b42cfa25 858
a86176e3 859 if(repa->GetUseChargeMatchingInClusterFinderSSD()==kTRUE) {
308b5ea4 860
861
a86176e3 862 //
863 // sign gold tracks
864 //
865 for (Int_t ip=0;ip<np;ip++){
8be4e1b1 866 Float_t xbest=1000,zbest=1000,qbest=0;
a86176e3 867 //
868 // select gold clusters
869 if ( (cnegative[ip]==1) && cpositive[negativepair[10*ip]]==1){
d036ccd3 870 Float_t yp=pos[ip].GetY();
a86176e3 871 Int_t j = negativepair[10*ip];
a64f9843 872
873 if( (pos[ip].GetQ()==0) && (neg[j].GetQ() ==0) ) {
874 // both bad, hence continue;
875 // mark both as used (to avoid recover at the end)
876 cused1[ip]++;
877 cused2[j]++;
878 continue;
879 }
880
a86176e3 881 ratio = (pos[ip].GetQ()-neg[j].GetQ())/(pos[ip].GetQ()+neg[j].GetQ());
d695268b 882 //cout<<"ratio="<<ratio<<endl;
a64f9843 883
884 // charge matching (note that if posQ or negQ is 0 -> ratio=1 and the following condition is met
2069484c 885 if (TMath::Abs(ratio)>0.2) continue; // note: 0.2=3xsigma_ratio calculated in cosmics tests
a64f9843 886
a86176e3 887 //
8be4e1b1 888 Float_t yn=neg[j].GetY();
889
890 Float_t xt, zt;
891 seg->GetPadCxz(yn, yp, xt, zt);
892
893 xbest=xt; zbest=zt;
2069484c 894
2069484c 895
a86176e3 896 qbest=0.5*(pos[ip].GetQ()+neg[j].GetQ());
d695268b 897 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 898
a86176e3 899 {
8be4e1b1 900 Double_t loc[3]={xbest,0.,zbest},trk[3]={0.,0.,0.};
a86176e3 901 mT2L->MasterToLocal(loc,trk);
902 lp[0]=trk[1];
903 lp[1]=trk[2];
00a7cc50 904 }
a86176e3 905 lp[4]=qbest; //Q
906 for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2;
907 for (Int_t ilab=0;ilab<3;ilab++){
908 milab[ilab] = pos[ip].GetLabel(ilab);
909 milab[ilab+3] = neg[j].GetLabel(ilab);
00a7cc50 910 }
04366a57 911 //
a86176e3 912 CheckLabels2(milab);
913 milab[3]=(((ip<<10) + j)<<10) + idet; // pos|neg|det
914 Int_t info[3] = {pos[ip].GetNd(),neg[j].GetNd(),fNlayer[fModule]};
d695268b 915
7eb157d7 916 lp[2]=4.968e-06; // 0.00223*0.00223; //SigmaY2
917 lp[3]=0.012; // 0.110*0.110; //SigmaZ2
d695268b 918 // out-of-diagonal element of covariance matrix
919 if( (info[0]==1) && (info[1]==1) ) lp[5]=-0.00012;
920 else if ( (info[0]>1) && (info[1]>1) ) {
7eb157d7 921 lp[2]=2.63e-06; // 0.0016*0.0016; //SigmaY2
922 lp[3]=0.0065; // 0.08*0.08; //SigmaZ2
923 lp[5]=-6.48e-05;
d695268b 924 }
925 else {
7eb157d7 926 lp[2]=4.80e-06; // 0.00219*0.00219
927 lp[3]=0.0093; // 0.0964*0.0964;
928 if (info[0]==1) {
929 lp[5]=-0.00014;
930 }
931 else {
932 lp[2]=2.79e-06; // 0.0017*0.0017;
933 lp[3]=0.00935; // 0.967*0.967;
934 lp[5]=-4.32e-05;
935 }
d695268b 936 }
937
a86176e3 938 AliITSRecPoint * cl2;
965c81a7 939 cl2 = new ((*clusters)[ncl]) AliITSRecPoint(milab,lp,info);
a86176e3 940
138df073 941 cl2->SetChargeRatio(ratio);
942 cl2->SetType(1);
943 fgPairs[ip*nn+j]=1;
a64f9843 944
138df073 945 if ((pos[ip].GetNd()+neg[j].GetNd())>6){ //multi cluster
946 cl2->SetType(2);
947 fgPairs[ip*nn+j]=2;
948 }
a64f9843 949
138df073 950 if(pos[ip].GetQ()==0) cl2->SetType(3);
951 if(neg[j].GetQ()==0) cl2->SetType(4);
a64f9843 952
138df073 953 cused1[ip]++;
954 cused2[j]++;
a64f9843 955
138df073 956 ncl++;
a86176e3 957 }
958 }
959
960 for (Int_t ip=0;ip<np;ip++){
8be4e1b1 961 Float_t xbest=1000,zbest=1000,qbest=0;
a86176e3 962 //
963 //
964 // select "silber" cluster
965 if ( cnegative[ip]==1 && cpositive[negativepair[10*ip]]==2){
966 Int_t in = negativepair[10*ip];
967 Int_t ip2 = positivepair[10*in];
968 if (ip2==ip) ip2 = positivepair[10*in+1];
969 Float_t pcharge = pos[ip].GetQ()+pos[ip2].GetQ();
a64f9843 970
d695268b 971
972
973 ratio = (pcharge-neg[in].GetQ())/(pcharge+neg[in].GetQ());
974 if ( (TMath::Abs(ratio)<0.2) && (pcharge!=0) ) {
975 //if ( (TMath::Abs(pcharge-neg[in].GetQ())<30) && (pcharge!=0) ) { //
a64f9843 976
a86176e3 977 //
978 // add first pair
a64f9843 979 if ( (fgPairs[ip*nn+in]==100)&&(pos[ip].GetQ() ) ) { //
980
d036ccd3 981 Float_t yp=pos[ip].GetY();
8be4e1b1 982 Float_t yn=neg[in].GetY();
983
984 Float_t xt, zt;
985 seg->GetPadCxz(yn, yp, xt, zt);
986
987 xbest=xt; zbest=zt;
2069484c 988
a86176e3 989 qbest =pos[ip].GetQ();
8be4e1b1 990 Double_t loc[3]={xbest,0.,zbest},trk[3]={0.,0.,0.};
a64f9843 991 mT2L->MasterToLocal(loc,trk);
992 lp[0]=trk[1];
993 lp[1]=trk[2];
994
a86176e3 995 lp[4]=qbest; //Q
996 for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2;
997 for (Int_t ilab=0;ilab<3;ilab++){
998 milab[ilab] = pos[ip].GetLabel(ilab);
999 milab[ilab+3] = neg[in].GetLabel(ilab);
1000 }
1001 //
1002 CheckLabels2(milab);
1003 ratio = (pos[ip].GetQ()-neg[in].GetQ())/(pos[ip].GetQ()+neg[in].GetQ());
1004 milab[3]=(((ip<<10) + in)<<10) + idet; // pos|neg|det
1005 Int_t info[3] = {pos[ip].GetNd(),neg[in].GetNd(),fNlayer[fModule]};
1006
7eb157d7 1007 lp[2]=4.968e-06; // 0.00223*0.00223; //SigmaY2
1008 lp[3]=0.012; // 0.110*0.110; //SigmaZ2
d695268b 1009 // out-of-diagonal element of covariance matrix
1010 if( (info[0]==1) && (info[1]==1) ) lp[5]=-0.00012;
1011 else if ( (info[0]>1) && (info[1]>1) ) {
7eb157d7 1012 lp[2]=2.63e-06; // 0.0016*0.0016; //SigmaY2
1013 lp[3]=0.0065; // 0.08*0.08; //SigmaZ2
1014 lp[5]=-6.48e-05;
d695268b 1015 }
1016 else {
7eb157d7 1017 lp[2]=4.80e-06; // 0.00219*0.00219
1018 lp[3]=0.0093; // 0.0964*0.0964;
1019 if (info[0]==1) {
1020 lp[5]=-0.00014;
1021 }
1022 else {
1023 lp[2]=2.79e-06; // 0.0017*0.0017;
1024 lp[3]=0.00935; // 0.967*0.967;
1025 lp[5]=-4.32e-05;
1026 }
d695268b 1027 }
1028
7eb157d7 1029 AliITSRecPoint * cl2;
138df073 1030 cl2 = new ((*clusters)[ncl]) AliITSRecPoint(milab,lp,info);
a86176e3 1031 cl2->SetChargeRatio(ratio);
1032 cl2->SetType(5);
1033 fgPairs[ip*nn+in] = 5;
1034 if ((pos[ip].GetNd()+neg[in].GetNd())>6){ //multi cluster
1035 cl2->SetType(6);
1036 fgPairs[ip*nn+in] = 6;
1037 }
a86176e3 1038 ncl++;
04366a57 1039 }
04366a57 1040
a64f9843 1041
04366a57 1042 //
a86176e3 1043 // add second pair
00a7cc50 1044
a86176e3 1045 // if (!(cused1[ip2] || cused2[in])){ //
a64f9843 1046 if ( (fgPairs[ip2*nn+in]==100) && (pos[ip2].GetQ()) ) {
1047
d036ccd3 1048 Float_t yp=pos[ip2].GetY();
8be4e1b1 1049 Float_t yn=neg[in].GetY();
1050
1051 Float_t xt, zt;
1052 seg->GetPadCxz(yn, yp, xt, zt);
1053
1054 xbest=xt; zbest=zt;
2069484c 1055
a86176e3 1056 qbest =pos[ip2].GetQ();
a64f9843 1057
8be4e1b1 1058 Double_t loc[3]={xbest,0.,zbest},trk[3]={0.,0.,0.};
a64f9843 1059 mT2L->MasterToLocal(loc,trk);
1060 lp[0]=trk[1];
1061 lp[1]=trk[2];
1062
a86176e3 1063 lp[4]=qbest; //Q
1064 for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2;
1065 for (Int_t ilab=0;ilab<3;ilab++){
1066 milab[ilab] = pos[ip2].GetLabel(ilab);
1067 milab[ilab+3] = neg[in].GetLabel(ilab);
00a7cc50 1068 }
a86176e3 1069 //
1070 CheckLabels2(milab);
1071 ratio = (pos[ip2].GetQ()-neg[in].GetQ())/(pos[ip2].GetQ()+neg[in].GetQ());
1072 milab[3]=(((ip2<<10) + in)<<10) + idet; // pos|neg|det
1073 Int_t info[3] = {pos[ip2].GetNd(),neg[in].GetNd(),fNlayer[fModule]};
308b5ea4 1074
7eb157d7 1075 lp[2]=4.968e-06; // 0.00223*0.00223; //SigmaY2
1076 lp[3]=0.012; // 0.110*0.110; //SigmaZ2
d695268b 1077 // out-of-diagonal element of covariance matrix
1078 if( (info[0]==1) && (info[1]==1) ) lp[5]=-0.00012;
1079 else if ( (info[0]>1) && (info[1]>1) ) {
7eb157d7 1080 lp[2]=2.63e-06; // 0.0016*0.0016; //SigmaY2
1081 lp[3]=0.0065; // 0.08*0.08; //SigmaZ2
1082 lp[5]=-6.48e-05;
d695268b 1083 }
1084 else {
7eb157d7 1085 lp[2]=4.80e-06; // 0.00219*0.00219
1086 lp[3]=0.0093; // 0.0964*0.0964;
1087 if (info[0]==1) {
1088 lp[5]=-0.00014;
1089 }
1090 else {
1091 lp[2]=2.79e-06; // 0.0017*0.0017;
1092 lp[3]=0.00935; // 0.967*0.967;
1093 lp[5]=-4.32e-05;
1094 }
d695268b 1095 }
1096
a86176e3 1097 AliITSRecPoint * cl2;
138df073 1098 cl2 = new ((*clusters)[ncl]) AliITSRecPoint(milab,lp,info);
a86176e3 1099
a86176e3 1100 cl2->SetChargeRatio(ratio);
1101 cl2->SetType(5);
1102 fgPairs[ip2*nn+in] =5;
1103 if ((pos[ip2].GetNd()+neg[in].GetNd())>6){ //multi cluster
1104 cl2->SetType(6);
1105 fgPairs[ip2*nn+in] =6;
1106 }
a86176e3 1107 ncl++;
a64f9843 1108 }
1109
a86176e3 1110 cused1[ip]++;
1111 cused1[ip2]++;
1112 cused2[in]++;
a64f9843 1113
1114 } // charge matching condition
1115
1116 } // 2 Pside cross 1 Nside
1117 } // loop over Pside clusters
a86176e3 1118
a64f9843 1119
1120
1121 //
1122 for (Int_t jn=0;jn<nn;jn++){
1123 if (cused2[jn]) continue;
8be4e1b1 1124 Float_t xbest=1000,zbest=1000,qbest=0;
a64f9843 1125 // select "silber" cluster
1126 if ( cpositive[jn]==1 && cnegative[positivepair[10*jn]]==2){
1127 Int_t ip = positivepair[10*jn];
1128 Int_t jn2 = negativepair[10*ip];
1129 if (jn2==jn) jn2 = negativepair[10*ip+1];
1130 Float_t pcharge = neg[jn].GetQ()+neg[jn2].GetQ();
04366a57 1131 //
a64f9843 1132
d695268b 1133
1134 ratio = (pcharge-pos[ip].GetQ())/(pcharge+pos[ip].GetQ());
1135 if ( (TMath::Abs(ratio)<0.2) && (pcharge!=0) ) {
1136
1137 /*
2069484c 1138 if ( (TMath::Abs(pcharge-pos[ip].GetQ())<30) && // charge matching
a64f9843 1139 (pcharge!=0) ) { // reject combinations of bad strips
d695268b 1140 */
1141
1142
a64f9843 1143 //
1144 // add first pair
1145 // if (!(cused1[ip]||cused2[jn])){
1146 if ( (fgPairs[ip*nn+jn]==100) && (neg[jn].GetQ()) ) { //
1147
d036ccd3 1148 Float_t yn=neg[jn].GetY();
8be4e1b1 1149 Float_t yp=pos[ip].GetY();
2069484c 1150
8be4e1b1 1151 Float_t xt, zt;
1152 seg->GetPadCxz(yn, yp, xt, zt);
1153
1154 xbest=xt; zbest=zt;
2069484c 1155
a64f9843 1156 qbest =neg[jn].GetQ();
d036ccd3 1157
a64f9843 1158 {
8be4e1b1 1159 Double_t loc[3]={xbest,0.,zbest},trk[3]={0.,0.,0.};
a64f9843 1160 mT2L->MasterToLocal(loc,trk);
1161 lp[0]=trk[1];
1162 lp[1]=trk[2];
b4704be3 1163 }
04366a57 1164
1165 lp[4]=qbest; //Q
1166 for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2;
1167 for (Int_t ilab=0;ilab<3;ilab++){
1168 milab[ilab] = pos[ip].GetLabel(ilab);
1169 milab[ilab+3] = neg[jn].GetLabel(ilab);
1170 }
1171 //
1172 CheckLabels2(milab);
1173 ratio = (pos[ip].GetQ()-neg[jn].GetQ())/(pos[ip].GetQ()+neg[jn].GetQ());
1174 milab[3]=(((ip<<10) + jn)<<10) + idet; // pos|neg|det
1175 Int_t info[3] = {pos[ip].GetNd(),neg[jn].GetNd(),fNlayer[fModule]};
1176
7eb157d7 1177 lp[2]=4.968e-06; // 0.00223*0.00223; //SigmaY2
1178 lp[3]=0.012; // 0.110*0.110; //SigmaZ2
d695268b 1179 // out-of-diagonal element of covariance matrix
1180 if( (info[0]==1) && (info[1]==1) ) lp[5]=-0.00012;
1181 else if ( (info[0]>1) && (info[1]>1) ) {
7eb157d7 1182 lp[2]=2.63e-06; // 0.0016*0.0016; //SigmaY2
1183 lp[3]=0.0065; // 0.08*0.08; //SigmaZ2
1184 lp[5]=-6.48e-05;
d695268b 1185 }
1186 else {
7eb157d7 1187 lp[2]=4.80e-06; // 0.00219*0.00219
1188 lp[3]=0.0093; // 0.0964*0.0964;
1189 if (info[0]==1) {
1190 lp[5]=-0.00014;
1191 }
1192 else {
1193 lp[2]=2.79e-06; // 0.0017*0.0017;
1194 lp[3]=0.00935; // 0.967*0.967;
1195 lp[5]=-4.32e-05;
1196 }
d695268b 1197 }
1198
00a7cc50 1199 AliITSRecPoint * cl2;
138df073 1200 cl2 = new ((*clusters)[ncl]) AliITSRecPoint(milab,lp,info);
308b5ea4 1201
00a7cc50 1202 cl2->SetChargeRatio(ratio);
1203 cl2->SetType(7);
308b5ea4 1204 fgPairs[ip*nn+jn] =7;
00a7cc50 1205 if ((pos[ip].GetNd()+neg[jn].GetNd())>6){ //multi cluster
1206 cl2->SetType(8);
308b5ea4 1207 fgPairs[ip*nn+jn]=8;
00a7cc50 1208 }
04366a57 1209 ncl++;
04366a57 1210 }
1211 //
1212 // add second pair
1213 // if (!(cused1[ip]||cused2[jn2])){
a64f9843 1214 if ( (fgPairs[ip*nn+jn2]==100)&&(neg[jn2].GetQ() ) ) { //
1215
d036ccd3 1216 Float_t yn=neg[jn2].GetY();
1217 Double_t yp=pos[ip].GetY();
2069484c 1218
8be4e1b1 1219 Float_t xt, zt;
1220 seg->GetPadCxz(yn, yp, xt, zt);
1221
1222 xbest=xt; zbest=zt;
2069484c 1223
04366a57 1224 qbest =neg[jn2].GetQ();
d036ccd3 1225
b4704be3 1226 {
8be4e1b1 1227 Double_t loc[3]={xbest,0.,zbest},trk[3]={0.,0.,0.};
b4704be3 1228 mT2L->MasterToLocal(loc,trk);
1229 lp[0]=trk[1];
1230 lp[1]=trk[2];
1231 }
d695268b 1232
04366a57 1233 lp[4]=qbest; //Q
1234 for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2;
1235 for (Int_t ilab=0;ilab<3;ilab++){
1236 milab[ilab] = pos[ip].GetLabel(ilab);
1237 milab[ilab+3] = neg[jn2].GetLabel(ilab);
1238 }
1239 //
1240 CheckLabels2(milab);
1241 ratio = (pos[ip].GetQ()-neg[jn2].GetQ())/(pos[ip].GetQ()+neg[jn2].GetQ());
1242 milab[3]=(((ip<<10) + jn2)<<10) + idet; // pos|neg|det
1243 Int_t info[3] = {pos[ip].GetNd(),neg[jn2].GetNd(),fNlayer[fModule]};
d695268b 1244
7eb157d7 1245 lp[2]=4.968e-06; // 0.00223*0.00223; //SigmaY2
1246 lp[3]=0.012; // 0.110*0.110; //SigmaZ2
d695268b 1247 // out-of-diagonal element of covariance matrix
1248 if( (info[0]==1) && (info[1]==1) ) lp[5]=-0.00012;
1249 else if ( (info[0]>1) && (info[1]>1) ) {
7eb157d7 1250 lp[2]=2.63e-06; // 0.0016*0.0016; //SigmaY2
1251 lp[3]=0.0065; // 0.08*0.08; //SigmaZ2
1252 lp[5]=-6.48e-05;
d695268b 1253 }
1254 else {
7eb157d7 1255 lp[2]=4.80e-06; // 0.00219*0.00219
1256 lp[3]=0.0093; // 0.0964*0.0964;
1257 if (info[0]==1) {
1258 lp[5]=-0.00014;
1259 }
1260 else {
1261 lp[2]=2.79e-06; // 0.0017*0.0017;
1262 lp[3]=0.00935; // 0.967*0.967;
1263 lp[5]=-4.32e-05;
1264 }
d695268b 1265 }
1266
00a7cc50 1267 AliITSRecPoint * cl2;
138df073 1268 cl2 = new ((*clusters)[ncl]) AliITSRecPoint(milab,lp,info);
308b5ea4 1269
308b5ea4 1270
00a7cc50 1271 cl2->SetChargeRatio(ratio);
308b5ea4 1272 fgPairs[ip*nn+jn2]=7;
00a7cc50 1273 cl2->SetType(7);
1274 if ((pos[ip].GetNd()+neg[jn2].GetNd())>6){ //multi cluster
1275 cl2->SetType(8);
308b5ea4 1276 fgPairs[ip*nn+jn2]=8;
00a7cc50 1277 }
04366a57 1278 ncl++;
04366a57 1279 }
1280 cused1[ip]++;
1281 cused2[jn]++;
1282 cused2[jn2]++;
308b5ea4 1283
a64f9843 1284 } // charge matching condition
1285
1286 } // 2 Nside cross 1 Pside
1287 } // loop over Pside clusters
1288
1289
1290
1291 for (Int_t ip=0;ip<np;ip++){
d695268b 1292
1293 if(cused1[ip]) continue;
1294
1295
8be4e1b1 1296 Float_t xbest=1000,zbest=1000,qbest=0;
04366a57 1297 //
a64f9843 1298 // 2x2 clusters
04366a57 1299 //
d695268b 1300 if ( (cnegative[ip]==2) && cpositive[negativepair[10*ip]]==2){
1301 Float_t minchargediff =4.;
15be6f41 1302 // Float_t minchargeratio =0.2;
d695268b 1303
1304 Int_t j=-1;
1305 for (Int_t di=0;di<cnegative[ip];di++){
1306 Int_t jc = negativepair[ip*10+di];
1307 Float_t chargedif = pos[ip].GetQ()-neg[jc].GetQ();
1308 ratio = (pos[ip].GetQ()-neg[jc].GetQ())/(pos[ip].GetQ()+neg[jc].GetQ());
1309 //if (TMath::Abs(chargedif)<minchargediff){
1310 if (TMath::Abs(ratio)<0.2){
1311 j =jc;
1312 minchargediff = TMath::Abs(chargedif);
15be6f41 1313 // minchargeratio = TMath::Abs(ratio);
d695268b 1314 }
1315 }
1316 if (j<0) continue; // not proper cluster
1317
1318
1319 Int_t count =0;
1320 for (Int_t di=0;di<cnegative[ip];di++){
1321 Int_t jc = negativepair[ip*10+di];
1322 Float_t chargedif = pos[ip].GetQ()-neg[jc].GetQ();
1323 if (TMath::Abs(chargedif)<minchargediff+3.) count++;
1324 }
1325 if (count>1) continue; // more than one "proper" cluster for positive
1326 //
1327
1328 count =0;
1329 for (Int_t dj=0;dj<cpositive[j];dj++){
1330 Int_t ic = positivepair[j*10+dj];
1331 Float_t chargedif = pos[ic].GetQ()-neg[j].GetQ();
1332 if (TMath::Abs(chargedif)<minchargediff+3.) count++;
1333 }
1334 if (count>1) continue; // more than one "proper" cluster for negative
1335
1336 Int_t jp = 0;
1337
1338 count =0;
1339 for (Int_t dj=0;dj<cnegative[jp];dj++){
1340 Int_t ic = positivepair[jp*10+dj];
1341 Float_t chargedif = pos[ic].GetQ()-neg[jp].GetQ();
1342 if (TMath::Abs(chargedif)<minchargediff+4.) count++;
1343 }
1344 if (count>1) continue;
1345 if (fgPairs[ip*nn+j]<100) continue;
1346 //
1347
1348
1349
1350 //almost gold clusters
1351 Float_t yp=pos[ip].GetY();
1352 Float_t yn=neg[j].GetY();
1353 Float_t xt, zt;
1354 seg->GetPadCxz(yn, yp, xt, zt);
1355 xbest=xt; zbest=zt;
1356 qbest=0.5*(pos[ip].GetQ()+neg[j].GetQ());
1357 {
1358 Double_t loc[3]={xbest,0.,zbest},trk[3]={0.,0.,0.};
1359 mT2L->MasterToLocal(loc,trk);
1360 lp[0]=trk[1];
1361 lp[1]=trk[2];
1362 }
1363 lp[4]=qbest; //Q
1364 for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2;
1365 for (Int_t ilab=0;ilab<3;ilab++){
1366 milab[ilab] = pos[ip].GetLabel(ilab);
1367 milab[ilab+3] = neg[j].GetLabel(ilab);
1368 }
1369 //
1370 CheckLabels2(milab);
1371 if ((neg[j].GetQ()==0)&&(pos[ip].GetQ()==0)) continue; // reject crosses of bad strips!!
1372 ratio = (pos[ip].GetQ()-neg[j].GetQ())/(pos[ip].GetQ()+neg[j].GetQ());
1373 milab[3]=(((ip<<10) + j)<<10) + idet; // pos|neg|det
1374 Int_t info[3] = {pos[ip].GetNd(),neg[j].GetNd(),fNlayer[fModule]};
1375
7eb157d7 1376 lp[2]=4.968e-06; // 0.00223*0.00223; //SigmaY2
1377 lp[3]=0.012; // 0.110*0.110; //SigmaZ2
d695268b 1378 // out-of-diagonal element of covariance matrix
1379 if( (info[0]==1) && (info[1]==1) ) lp[5]=-0.00012;
1380 else if ( (info[0]>1) && (info[1]>1) ) {
7eb157d7 1381 lp[2]=2.63e-06; // 0.0016*0.0016; //SigmaY2
1382 lp[3]=0.0065; // 0.08*0.08; //SigmaZ2
1383 lp[5]=-6.48e-05;
d695268b 1384 }
1385 else {
7eb157d7 1386 lp[2]=4.80e-06; // 0.00219*0.00219
1387 lp[3]=0.0093; // 0.0964*0.0964;
1388 if (info[0]==1) {
1389 lp[5]=-0.00014;
1390 }
1391 else {
1392 lp[2]=2.79e-06; // 0.0017*0.0017;
1393 lp[3]=0.00935; // 0.967*0.967;
1394 lp[5]=-4.32e-05;
1395 }
d695268b 1396 }
1397
1398 AliITSRecPoint * cl2;
138df073 1399 cl2 = new ((*clusters)[ncl]) AliITSRecPoint(milab,lp,info);
d695268b 1400
1401 cl2->SetChargeRatio(ratio);
1402 cl2->SetType(10);
1403 fgPairs[ip*nn+j]=10;
1404 if ((pos[ip].GetNd()+neg[j].GetNd())>6){ //multi cluster
1405 cl2->SetType(11);
1406 fgPairs[ip*nn+j]=11;
1407 }
1408 cused1[ip]++;
1409 cused2[j]++;
d695268b 1410 ncl++;
1411
1412 } // 2X2
1413 } // loop over Pside 1Dclusters
1414
1415
1416
1417 for (Int_t ip=0;ip<np;ip++){
1418
1419 if(cused1[ip]) continue;
1420
1421
1422 Float_t xbest=1000,zbest=1000,qbest=0;
1423 //
1424 // manyxmany clusters
1425 //
a64f9843 1426 if ( (cnegative[ip]<5) && cpositive[negativepair[10*ip]]<5){
1427 Float_t minchargediff =4.;
1428 Int_t j=-1;
1429 for (Int_t di=0;di<cnegative[ip];di++){
1430 Int_t jc = negativepair[ip*10+di];
1431 Float_t chargedif = pos[ip].GetQ()-neg[jc].GetQ();
1432 if (TMath::Abs(chargedif)<minchargediff){
1433 j =jc;
1434 minchargediff = TMath::Abs(chargedif);
1435 }
00a7cc50 1436 }
a64f9843 1437 if (j<0) continue; // not proper cluster
1438
1439 Int_t count =0;
1440 for (Int_t di=0;di<cnegative[ip];di++){
1441 Int_t jc = negativepair[ip*10+di];
1442 Float_t chargedif = pos[ip].GetQ()-neg[jc].GetQ();
1443 if (TMath::Abs(chargedif)<minchargediff+3.) count++;
00a7cc50 1444 }
a64f9843 1445 if (count>1) continue; // more than one "proper" cluster for positive
1446 //
00a7cc50 1447
a64f9843 1448 count =0;
1449 for (Int_t dj=0;dj<cpositive[j];dj++){
1450 Int_t ic = positivepair[j*10+dj];
1451 Float_t chargedif = pos[ic].GetQ()-neg[j].GetQ();
1452 if (TMath::Abs(chargedif)<minchargediff+3.) count++;
1453 }
1454 if (count>1) continue; // more than one "proper" cluster for negative
1455
1456 Int_t jp = 0;
1457
1458 count =0;
1459 for (Int_t dj=0;dj<cnegative[jp];dj++){
1460 Int_t ic = positivepair[jp*10+dj];
1461 Float_t chargedif = pos[ic].GetQ()-neg[jp].GetQ();
1462 if (TMath::Abs(chargedif)<minchargediff+4.) count++;
1463 }
1464 if (count>1) continue;
1465 if (fgPairs[ip*nn+j]<100) continue;
1466 //
1467
1468 //almost gold clusters
d036ccd3 1469 Float_t yp=pos[ip].GetY();
1470 Float_t yn=neg[j].GetY();
1471
8be4e1b1 1472
1473 Float_t xt, zt;
1474 seg->GetPadCxz(yn, yp, xt, zt);
1475
1476 xbest=xt; zbest=zt;
1477
a64f9843 1478 qbest=0.5*(pos[ip].GetQ()+neg[j].GetQ());
d036ccd3 1479
a64f9843 1480 {
8be4e1b1 1481 Double_t loc[3]={xbest,0.,zbest},trk[3]={0.,0.,0.};
a64f9843 1482 mT2L->MasterToLocal(loc,trk);
1483 lp[0]=trk[1];
1484 lp[1]=trk[2];
1485 }
a64f9843 1486 lp[4]=qbest; //Q
1487 for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2;
1488 for (Int_t ilab=0;ilab<3;ilab++){
1489 milab[ilab] = pos[ip].GetLabel(ilab);
1490 milab[ilab+3] = neg[j].GetLabel(ilab);
1491 }
1492 //
1493 CheckLabels2(milab);
db6e54cd 1494 if ((neg[j].GetQ()==0)&&(pos[ip].GetQ()==0)) continue; // reject crosses of bad strips!!
a64f9843 1495 ratio = (pos[ip].GetQ()-neg[j].GetQ())/(pos[ip].GetQ()+neg[j].GetQ());
1496 milab[3]=(((ip<<10) + j)<<10) + idet; // pos|neg|det
1497 Int_t info[3] = {pos[ip].GetNd(),neg[j].GetNd(),fNlayer[fModule]};
d695268b 1498
7eb157d7 1499 lp[2]=4.968e-06; // 0.00223*0.00223; //SigmaY2
1500 lp[3]=0.012; // 0.110*0.110; //SigmaZ2
d695268b 1501 // out-of-diagonal element of covariance matrix
1502 if( (info[0]==1) && (info[1]==1) ) lp[5]=-0.00012;
1503 else if ( (info[0]>1) && (info[1]>1) ) {
7eb157d7 1504 lp[2]=2.63e-06; // 0.0016*0.0016; //SigmaY2
1505 lp[3]=0.0065; // 0.08*0.08; //SigmaZ2
1506 lp[5]=-6.48e-05;
d695268b 1507 }
1508 else {
7eb157d7 1509 lp[2]=4.80e-06; // 0.00219*0.00219
1510 lp[3]=0.0093; // 0.0964*0.0964;
1511 if (info[0]==1) {
1512 lp[5]=-0.00014;
1513 }
1514 else {
1515 lp[2]=2.79e-06; // 0.0017*0.0017;
1516 lp[3]=0.00935; // 0.967*0.967;
1517 lp[5]=-4.32e-05;
1518 }
d695268b 1519 }
6c8e94cf 1520 //
1521 if (fRawID2ClusID) { // set rawID <-> clusterID correspondence for embedding
1522 const int kMaxRefRW = 200;
1523 UInt_t nrefsRW,refsRW[kMaxRefRW];
1524 nrefsRW = fRawIDRef[0].GetReferences(j,refsRW,kMaxRefRW); // n-side
1525 for (int ir=nrefsRW;ir--;) {
1526 int rwid = (int)refsRW[ir];
1527 if (fRawID2ClusID->GetSize()<=rwid) fRawID2ClusID->Set( (rwid+10)<<1 );
1528 (*fRawID2ClusID)[rwid] = fNClusters+1; // RS: store clID+1 as a reference to the cluster
1529 }
1530 //
1531 nrefsRW = fRawIDRef[1].GetReferences(ip,refsRW,kMaxRefRW); // p-side
1532 for (int ir=nrefsRW;ir--;) {
1533 int rwid = (int)refsRW[ir];
1534 if (fRawID2ClusID->GetSize()<=rwid) fRawID2ClusID->Set( (rwid+10)<<1 );
1535 (*fRawID2ClusID)[rwid] = fNClusters+1; // RS: store clID+1 as a reference to the cluster
1536 }
1537 //
1538 milab[0] = fNClusters+1; // RS: assign id as cluster label
1539 }
d695268b 1540
a64f9843 1541 AliITSRecPoint * cl2;
138df073 1542 cl2 = new ((*clusters)[ncl]) AliITSRecPoint(milab,lp,info);
d036ccd3 1543
a64f9843 1544 cl2->SetChargeRatio(ratio);
d695268b 1545 cl2->SetType(12);
1546 fgPairs[ip*nn+j]=12;
a64f9843 1547 if ((pos[ip].GetNd()+neg[j].GetNd())>6){ //multi cluster
d695268b 1548 cl2->SetType(13);
1549 fgPairs[ip*nn+j]=13;
a64f9843 1550 }
1551 cused1[ip]++;
1552 cused2[j]++;
6c8e94cf 1553 ncl++;
1554 fNClusters++;
a64f9843 1555
1556 } // manyXmany
1557 } // loop over Pside 1Dclusters
1558
a64f9843 1559 } // use charge matching
1560
04366a57 1561
a64f9843 1562 // recover all the other crosses
04366a57 1563 //
1564 for (Int_t i=0; i<np; i++) {
8be4e1b1 1565 Float_t xbest=1000,zbest=1000,qbest=0;
d036ccd3 1566 Float_t yp=pos[i].GetY();
a64f9843 1567 if ((pos[i].GetQ()>0)&&(pos[i].GetQ()<3)) continue;
04366a57 1568 for (Int_t j=0; j<nn; j++) {
1569 // for (Int_t di = 0;di<cpositive[i];di++){
1570 // Int_t j = negativepair[10*i+di];
a64f9843 1571 if ((neg[j].GetQ()>0)&&(neg[j].GetQ()<3)) continue;
1572
1573 if ((neg[j].GetQ()==0)&&(pos[i].GetQ()==0)) continue; // reject crosses of bad strips!!
1574
04366a57 1575 if (cused2[j]||cused1[i]) continue;
308b5ea4 1576 if (fgPairs[i*nn+j]>0 &&fgPairs[i*nn+j]<100) continue;
04366a57 1577 ratio = (pos[i].GetQ()-neg[j].GetQ())/(pos[i].GetQ()+neg[j].GetQ());
d036ccd3 1578 Float_t yn=neg[j].GetY();
1579
8be4e1b1 1580 Float_t xt, zt;
1581 seg->GetPadCxz(yn, yp, xt, zt);
d036ccd3 1582
5846520e 1583 if (TMath::Abs(xt)<hwSSD)
1584 if (TMath::Abs(zt)<hlSSD) {
8be4e1b1 1585 xbest=xt; zbest=zt;
1586
04366a57 1587 qbest=0.5*(pos[i].GetQ()+neg[j].GetQ());
d036ccd3 1588
b4704be3 1589 {
8be4e1b1 1590 Double_t loc[3]={xbest,0.,zbest},trk[3]={0.,0.,0.};
b4704be3 1591 mT2L->MasterToLocal(loc,trk);
1592 lp[0]=trk[1];
1593 lp[1]=trk[2];
1594 }
04366a57 1595 lp[4]=qbest; //Q
1596 for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2;
1597 for (Int_t ilab=0;ilab<3;ilab++){
1598 milab[ilab] = pos[i].GetLabel(ilab);
1599 milab[ilab+3] = neg[j].GetLabel(ilab);
1600 }
1601 //
1602 CheckLabels2(milab);
1603 milab[3]=(((i<<10) + j)<<10) + idet; // pos|neg|det
1604 Int_t info[3] = {pos[i].GetNd(),neg[j].GetNd(),fNlayer[fModule]};
d695268b 1605
7eb157d7 1606 lp[2]=4.968e-06; // 0.00223*0.00223; //SigmaY2
1607 lp[3]=0.012; // 0.110*0.110; //SigmaZ2
d695268b 1608 // out-of-diagonal element of covariance matrix
1609 if( (info[0]==1) && (info[1]==1) ) lp[5]=-0.00012;
1610 else if ( (info[0]>1) && (info[1]>1) ) {
7eb157d7 1611 lp[2]=2.63e-06; // 0.0016*0.0016; //SigmaY2
1612 lp[3]=0.0065; // 0.08*0.08; //SigmaZ2
1613 lp[5]=-6.48e-05;
d695268b 1614 }
1615 else {
7eb157d7 1616 lp[2]=4.80e-06; // 0.00219*0.00219
1617 lp[3]=0.0093; // 0.0964*0.0964;
1618 if (info[0]==1) {
1619 lp[5]=-0.00014;
1620 }
1621 else {
1622 lp[2]=2.79e-06; // 0.0017*0.0017;
1623 lp[3]=0.00935; // 0.967*0.967;
1624 lp[5]=-4.32e-05;
1625 }
d695268b 1626 }
1627
00a7cc50 1628 AliITSRecPoint * cl2;
138df073 1629 cl2 = new ((*clusters)[ncl]) AliITSRecPoint(milab,lp,info);
308b5ea4 1630
00a7cc50 1631 cl2->SetChargeRatio(ratio);
1632 cl2->SetType(100+cpositive[j]+cnegative[i]);
a64f9843 1633
1634 if(pos[i].GetQ()==0) cl2->SetType(200+cpositive[j]+cnegative[i]);
1635 if(neg[j].GetQ()==0) cl2->SetType(300+cpositive[j]+cnegative[i]);
04366a57 1636 ncl++;
04366a57 1637 }
1638 }
1639 }
d695268b 1640
1641
d695268b 1642
7101948c 1643 if(repa->GetUseBadChannelsInClusterFinderSSD()==kTRUE) {
1644
1645 //---------------------------------------------------------
1646 // recover crosses of good 1D clusters with bad strips on the other side
1647 // Note1: at first iteration skip modules with a bad side (or almost), (would produce too many fake!)
1648 // Note2: for modules with a bad side see below
1649
1650 AliITSCalibrationSSD* cal = (AliITSCalibrationSSD*)GetResp(fModule);
1651 Int_t countPbad=0, countNbad=0;
1652 for(Int_t ib=0; ib<768; ib++) {
1653 if(cal->IsPChannelBad(ib)) countPbad++;
1654 if(cal->IsNChannelBad(ib)) countNbad++;
1655 }
1656 // AliInfo(Form("module %d has %d P- and %d N-bad strips",fModule,countPbad,countNbad));
d695268b 1657
7101948c 1658 if( (countPbad<100) && (countNbad<100) ) { // no bad side!!
d695268b 1659
7101948c 1660 for (Int_t i=0; i<np; i++) { // loop over Nside 1Dclusters with no crosses
1661 if(cnegative[i]) continue; // if intersecting Pside clusters continue;
d695268b 1662
7101948c 1663 // for(Int_t ib=0; ib<768; ib++) { // loop over all Pstrips
1664 for(Int_t ib=15; ib<753; ib++) { // loop over all Pstrips
d695268b 1665
7101948c 1666 if(cal->IsPChannelBad(ib)) { // check if strips is bad
1667 Float_t yN=pos[i].GetY();
1668 Float_t xt, zt;
1669 seg->GetPadCxz(1.*ib, yN, xt, zt);
d695268b 1670
7101948c 1671 //----------
1672 // bad Pstrip is crossing the Nside 1Dcluster -> create recpoint
1673 //
5846520e 1674 if ( (TMath::Abs(xt)<hwSSD) && (TMath::Abs(zt)<hlSSD) ) {
7101948c 1675 Double_t loc[3]={xt,0.,zt},trk[3]={0.,0.,0.};
1676 mT2L->MasterToLocal(loc,trk);
1677 lp[0]=trk[1];
1678 lp[1]=trk[2];
1679 lp[4]=pos[i].GetQ(); //Q
1680 for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2;
1681 for (Int_t ilab=0;ilab<3;ilab++) milab[ilab] = pos[i].GetLabel(ilab);
1682 CheckLabels2(milab);
1683 milab[3]=( (i<<10) << 10 ) + idet; // pos|neg|det
1684 Int_t info[3] = {pos[i].GetNd(),0,fNlayer[fModule]};
1685
7eb157d7 1686 lp[2]=4.968e-06; // 0.00223*0.00223; //SigmaY2
1687 lp[3]=0.012; // 0.110*0.110; //SigmaZ2
7101948c 1688 lp[5]=-0.00012; // out-of-diagonal element of covariance matrix
7eb157d7 1689 if (info[0]>1) {
1690 lp[2]=4.80e-06;
1691 lp[3]=0.0093;
1692 lp[5]=0.00014;
1693 }
1694
7101948c 1695 AliITSRecPoint * cl2;
138df073 1696 cl2 = new ((*clusters)[ncl]) AliITSRecPoint(milab,lp,info);
7101948c 1697 cl2->SetChargeRatio(1.);
d695268b 1698 cl2->SetType(50);
7101948c 1699 ncl++;
1700 } // cross is within the detector
1701 //
1702 //--------------
1703
1704 } // bad Pstrip
d695268b 1705
7101948c 1706 } // end loop over Pstrips
d695268b 1707
7101948c 1708 } // end loop over Nside 1D clusters
d695268b 1709
7101948c 1710 for (Int_t j=0; j<nn; j++) { // loop over Pside 1D clusters with no crosses
1711 if(cpositive[j]) continue;
1712
1713 // for(Int_t ib=0; ib<768; ib++) { // loop over all Nside strips
1714 for(Int_t ib=15; ib<753; ib++) { // loop over all Nside strips
1715
1716 if(cal->IsNChannelBad(ib)) { // check if strip is bad
1717 Float_t yP=neg[j].GetY();
1718 Float_t xt, zt;
1719 seg->GetPadCxz(yP, 1.*ib, xt, zt);
1720
1721 //----------
1722 // bad Nstrip is crossing the Pside 1Dcluster -> create recpoint
1723 //
5846520e 1724 if ( (TMath::Abs(xt)<hwSSD) && (TMath::Abs(zt)<hlSSD) ) {
7101948c 1725 Double_t loc[3]={xt,0.,zt},trk[3]={0.,0.,0.};
1726 mT2L->MasterToLocal(loc,trk);
1727 lp[0]=trk[1];
1728 lp[1]=trk[2];
1729 lp[4]=neg[j].GetQ(); //Q
1730 for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2;
1731 for (Int_t ilab=0;ilab<3;ilab++) milab[ilab] = neg[j].GetLabel(ilab);
1732 CheckLabels2(milab);
1733 milab[3]=( j << 10 ) + idet; // pos|neg|det
1734 Int_t info[3]={0,(Int_t)neg[j].GetNd(),fNlayer[fModule]};
1735
7eb157d7 1736 lp[2]=4.968e-06; // 0.00223*0.00223; //SigmaY2
1737 lp[3]=0.012; // 0.110*0.110; //SigmaZ2
7101948c 1738 lp[5]=-0.00012; // out-of-diagonal element of covariance matrix
7eb157d7 1739 if (info[0]>1) {
1740 lp[2]=2.79e-06;
1741 lp[3]=0.00935;
1742 lp[5]=-4.32e-05;
1743 }
7101948c 1744
1745 AliITSRecPoint * cl2;
138df073 1746 cl2 = new ((*clusters)[ncl]) AliITSRecPoint(milab,lp,info);
7101948c 1747 cl2->SetChargeRatio(1.);
1748 cl2->SetType(60);
7101948c 1749 ncl++;
1750 } // cross is within the detector
1751 //
1752 //--------------
1753
1754 } // bad Nstrip
1755 } // end loop over Nstrips
1756 } // end loop over Pside 1D clusters
1757
1758 } // no bad sides
1759
1760 //---------------------------------------------------------
d695268b 1761
7101948c 1762 else if( (countPbad>700) && (countNbad<100) ) { // bad Pside!!
d695268b 1763
7101948c 1764 for (Int_t i=0; i<np; i++) { // loop over Nside 1Dclusters with no crosses
1765 if(cnegative[i]) continue; // if intersecting Pside clusters continue;
d695268b 1766
7101948c 1767 Float_t xt, zt;
1768 Float_t yN=pos[i].GetY();
1769 Float_t yP=0.;
1770 if (seg->GetLayer()==5) yP = yN + (7.6/1.9);
1771 else yP = yN - (7.6/1.9);
1772 seg->GetPadCxz(yP, yN, xt, zt);
1773
5846520e 1774 if ( (TMath::Abs(xt)<hwSSD) && (TMath::Abs(zt)<hlSSD) ) {
7101948c 1775 Double_t loc[3]={xt,0.,zt},trk[3]={0.,0.,0.};
1776 mT2L->MasterToLocal(loc,trk);
1777 lp[0]=trk[1];
1778 lp[1]=trk[2];
1779 lp[4]=pos[i].GetQ(); //Q
1780 for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2;
1781 for (Int_t ilab=0;ilab<3;ilab++) milab[ilab] = pos[i].GetLabel(ilab);
1782 CheckLabels2(milab);
1783 milab[3]=( (i<<10) << 10 ) + idet; // pos|neg|det
1784 Int_t info[3] = {(Int_t)pos[i].GetNd(),0,fNlayer[fModule]};
d695268b 1785
7eb157d7 1786 lp[2]=0.00098; // 0.031*0.031; //SigmaY2
1787 lp[3]=1.329; // 1.15*1.15; //SigmaZ2
1788 lp[5]=-0.0359;
1789 if(info[0]>1) lp[2]=0.00097;
1790
7101948c 1791 AliITSRecPoint * cl2;
138df073 1792 cl2 = new ((*clusters)[ncl]) AliITSRecPoint(milab,lp,info);
7101948c 1793 cl2->SetChargeRatio(1.);
1794 cl2->SetType(70);
7101948c 1795 ncl++;
1796 } // cross is within the detector
1797 //
1798 //--------------
1799
1800 } // end loop over Nside 1D clusters
1801
1802 } // bad Pside module
d695268b 1803
7101948c 1804 else if( (countNbad>700) && (countPbad<100) ) { // bad Nside!!
d695268b 1805
7101948c 1806 for (Int_t j=0; j<nn; j++) { // loop over Pside 1D clusters with no crosses
1807 if(cpositive[j]) continue;
1808
1809 Float_t xt, zt;
1810 Float_t yP=neg[j].GetY();
1811 Float_t yN=0.;
1812 if (seg->GetLayer()==5) yN = yP - (7.6/1.9);
1813 else yN = yP + (7.6/1.9);
1814 seg->GetPadCxz(yP, yN, xt, zt);
1815
5846520e 1816 if ( (TMath::Abs(xt)<hwSSD) && (TMath::Abs(zt)<hlSSD) ) {
7101948c 1817 Double_t loc[3]={xt,0.,zt},trk[3]={0.,0.,0.};
1818 mT2L->MasterToLocal(loc,trk);
1819 lp[0]=trk[1];
1820 lp[1]=trk[2];
1821 lp[4]=neg[j].GetQ(); //Q
1822 for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2;
1823 for (Int_t ilab=0;ilab<3;ilab++) milab[ilab] = neg[j].GetLabel(ilab);
1824 CheckLabels2(milab);
1825 milab[3]=( j << 10 ) + idet; // pos|neg|det
1826 Int_t info[3] = {0,(Int_t)neg[j].GetNd(),fNlayer[fModule]};
d695268b 1827
7eb157d7 1828 lp[2]=7.27e-05; // 0.0085*0.0085; //SigmaY2
1829 lp[3]=1.33; // 1.15*1.15; //SigmaZ2
1830 lp[5]=0.00931;
1831 if(info[1]>1) lp[2]=6.91e-05;
7101948c 1832
1833 AliITSRecPoint * cl2;
138df073 1834 cl2 = new ((*clusters)[ncl]) AliITSRecPoint(milab,lp,info);
7101948c 1835 cl2->SetChargeRatio(1.);
1836 cl2->SetType(80);
7101948c 1837 ncl++;
1838 } // cross is within the detector
1839 //
1840 //--------------
1841
1842 } // end loop over Pside 1D clusters
d695268b 1843
7101948c 1844 } // bad Nside module
d695268b 1845
7101948c 1846 //---------------------------------------------------------
d695268b 1847
7101948c 1848 } // use bad channels
d695268b 1849
d695268b 1850 //cout<<ncl<<" clusters for this module"<<endl;
1851
8be4e1b1 1852 delete [] cnegative;
1853 delete [] cused1;
1854 delete [] negativepair;
1855 delete [] cpositive;
1856 delete [] cused2;
1857 delete [] positivepair;
04366a57 1858
0a56760a 1859}