Line data Source code
1 : /**************************************************************************
2 : * Copyright(c) 2007-2009, ALICE Experiment at CERN, All rights reserved. *
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 : **************************************************************************/
15 :
16 : /* $Id: AliHLTITSClusterFinderSSD.cxx 34920 2009-09-22 07:48:53Z masera $ */
17 :
18 : ////////////////////////////////////////////////////////////////////////////
19 : // Implementation of the ITS clusterer V2 class //
20 : // //
21 : // Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch //
22 : // Last revision: 13-05-09 Enrico Fragiacomo //
23 : // enrico.fragiacomo@ts.infn.it //
24 : // //
25 : ///////////////////////////////////////////////////////////////////////////
26 :
27 : #include <Riostream.h>
28 : #include "AliLog.h"
29 :
30 : #include "AliHLTITSClusterFinderSSD.h"
31 : #include "AliITSRecPoint.h"
32 : #include "AliITSgeomTGeo.h"
33 : #include "AliITSDetTypeRec.h"
34 : #include "AliRawReader.h"
35 : #include "AliITSRawStreamSSD.h"
36 : #include <TClonesArray.h>
37 : #include "AliITSdigitSSD.h"
38 : #include "AliITSReconstructor.h"
39 : #include "AliITSCalibrationSSD.h"
40 : #include "AliITSsegmentationSSD.h"
41 :
42 : Short_t *AliHLTITSClusterFinderSSD::fgPairs = 0x0;
43 : Int_t AliHLTITSClusterFinderSSD::fgPairsSize = 0;
44 : const Float_t AliHLTITSClusterFinderSSD::fgkThreshold = 5.;
45 :
46 : const Float_t AliHLTITSClusterFinderSSD::fgkCosmic2008StripShifts[16][9] =
47 : {{-0.35,-0.35,-0.35,-0.35,-0.35,-0.35,-0.35,-0.35,-0.35}, // DDL 512
48 : {-0.35,-0.35,-0.35,-0.35,-0.35,-0.35,-0.35,-0.35,-0.35}, // DDL 513
49 : {-0.15,-0.15,-0.15,-0.15,-0.15,-0.15,-0.15,-0.15,-0.15}, // DDL 514
50 : {-0.15,-0.15,-0.15,-0.15,-0.15,-0.15,-0.15,-0.15,-0.15}, // DDL 515
51 : { 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00}, // DDL 516
52 : { 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00}, // DDL 517
53 : {-0.15,-0.15,-0.15,-0.15,-0.15,-0.15,-0.15,-0.15,-0.15}, // DDL 518
54 : {-0.15,-0.15,-0.15,-0.15,-0.15,-0.15,-0.15,-0.15,-0.15}, // DDL 519
55 : {-0.15,-0.15,-0.15,-0.15,-0.15,-0.15,-0.15,-0.25,-0.15}, // DDL 520
56 : {-0.15,-0.15,-0.15,-0.15,-0.15,-0.15,-0.15,-0.15,-0.15}, // DDL 521
57 : {-0.10,-0.10,-0.10,-0.40,-0.40,-0.40,-0.10,-0.10,-0.45}, // DDL 522
58 : {-0.10,-0.10,-0.10,-0.35,-0.35,-0.35,-0.10,-0.35,-0.50}, // DDL 523
59 : { 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00}, // DDL 524
60 : { 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00}, // DDL 525
61 : { 0.35, 0.35, 0.35, 0.35, 0.35, 0.35, 0.35, 0.35, 0.35}, // DDL 526
62 : { 0.45, 0.45, 0.45, 0.45, 0.45, 0.45, 0.45, 0.45, 0.45}}; // DDL 527
63 :
64 6 : ClassImp(AliHLTITSClusterFinderSSD)
65 :
66 :
67 : AliHLTITSClusterFinderSSD::AliHLTITSClusterFinderSSD(AliITSDetTypeRec* dettyp, AliRawReader *reader)
68 : :
69 0 : AliITSClusterFinder(dettyp),
70 0 : fRecoParam(0),
71 0 : fRawReader( reader ),
72 0 : fRawStream( 0 ),
73 0 : fLastSSD1(AliITSgeomTGeo::GetModuleIndex(6,1,1)-1)
74 0 : {
75 : //Default constructor
76 : //
77 : // Initialisation of ITS geometry
78 : //
79 0 : Int_t mmax=AliITSgeomTGeo::GetNModules();
80 0 : for (Int_t m=0; m<mmax; m++) {
81 0 : Int_t lay,lad,det; AliITSgeomTGeo::GetModuleId(m,lay,lad,det);
82 0 : fNdet[m] = (lad-1)*AliITSgeomTGeo::GetNDetectors(lay) + (det-1);
83 0 : fNlayer[m] = lay-1;
84 0 : }
85 :
86 0 : fRecoParam = (AliITSRecoParam*) AliITSReconstructor::GetRecoParam();
87 0 : if( !fRecoParam ){
88 0 : fRecoParam = AliITSRecoParam::GetHighFluxParam();
89 : // AliWarning("Using default AliITSRecoParam class");
90 0 : }
91 0 : fRawStream = new AliITSRawStreamSSD( fRawReader);
92 0 : }
93 :
94 : //______________________________________________________________________
95 0 : AliHLTITSClusterFinderSSD::AliHLTITSClusterFinderSSD(const AliHLTITSClusterFinderSSD &cf) : AliITSClusterFinder(cf), fRecoParam(cf.fRecoParam), fRawReader(cf.fRawReader), fRawStream(0), fLastSSD1(cf.fLastSSD1)
96 0 : {
97 : // Dummy
98 0 : }
99 :
100 : //______________________________________________________________________
101 : AliHLTITSClusterFinderSSD& AliHLTITSClusterFinderSSD::operator=(const AliHLTITSClusterFinderSSD& ){
102 : // Dummy
103 0 : return *this;
104 : }
105 :
106 : AliHLTITSClusterFinderSSD::~AliHLTITSClusterFinderSSD()
107 0 : {
108 : // destructor
109 0 : delete fRawStream;
110 0 : }
111 :
112 : void AliHLTITSClusterFinderSSD::RawdataToClusters( std::vector<AliITSRecPoint> &clusters )
113 : {
114 : //------------------------------------------------------------
115 : // Actual SSD cluster finder for raw data
116 : //------------------------------------------------------------
117 :
118 0 : fRawReader->Reset();
119 :
120 : const Int_t kNADC = 12;
121 : const Int_t kMaxADCClusters = 1000;
122 :
123 0 : Int_t strips[kNADC][2][kMaxADCClusters][2]; // [ADC],[side],[istrip], [0]=istrip [1]=signal
124 0 : Int_t nStrips[kNADC][2];
125 :
126 0 : for( int i=0; i<kNADC; i++ ){
127 0 : nStrips[i][0] = 0;
128 0 : nStrips[i][1] = 0;
129 : }
130 :
131 : Int_t ddl = -1;
132 : Int_t ad = -1;
133 :
134 : //*
135 : //* Loop over modules DDL+AD
136 : //*
137 :
138 0 : while (kTRUE) {
139 :
140 0 : bool next = fRawStream->Next();
141 :
142 : //*
143 : //* Continue if corrupted input
144 : //*
145 :
146 0 : if( (!next)&&(fRawStream->flag) ){
147 0 : AliWarning(Form("HLT ClustersFinderSSD: Corrupted data: warning from RawReader"));
148 0 : continue;
149 : }
150 :
151 0 : Int_t newDDL = fRawStream->GetDDL();
152 0 : Int_t newAD = fRawStream->GetAD();
153 :
154 0 : if( next ){
155 0 : if( newDDL<0 || newDDL>15 ){
156 : // AliWarning(Form("HLT ClustersFinderSSD: Corrupted data: wrong DDL number (%d)",newDDL));
157 0 : continue;
158 : }
159 :
160 0 : if( newAD<1 || newAD>9 ){
161 : // AliWarning(Form("HLT ClustersFinderSSD: Corrupted data: wrong AD number (%d)",newAD));
162 0 : continue;
163 : }
164 : }
165 :
166 0 : bool newModule = ( !next || ddl!= newDDL || ad!=newAD );
167 :
168 0 : if( newModule && ddl>=0 && ad>=0 ){
169 :
170 : //*
171 : //* Reconstruct the previous module --- actual clusterfinder
172 : //*
173 : //cout<<endl;
174 0 : for( int adc = 0; adc<kNADC; adc++ ){
175 :
176 : //* 1D clusterfinder
177 :
178 0 : Ali1Dcluster clusters1D[2][kMaxADCClusters]; // per ADC, per side
179 0 : Int_t nClusters1D[2] = {0,0};
180 : //int nstat[2] = {0,0};
181 0 : fModule = AliITSRawStreamSSD::GetModuleNumber(ddl, (ad - 1) * 12 + adc );
182 :
183 0 : if( fModule<0 ){
184 : // AliWarning(Form("HLT ClustersFinderSSD: Corrupted data: module (ddl %d ad %d adc %d) not found in the map",ddl,ad,adc));
185 0 : continue;
186 : }
187 :
188 0 : AliITSCalibrationSSD* cal = (AliITSCalibrationSSD*)fDetTypeRec->GetCalibrationModel(fModule);
189 0 : if( !cal ){
190 0 : AliWarning(Form("HLT ClustersFinderSSD: No calibration found for module (ddl %d ad %d adc %d)",ddl,ad,adc));
191 0 : continue;
192 : }
193 :
194 : Float_t dStrip = 0;
195 :
196 0 : if( fRecoParam->GetUseCosmicRunShiftsSSD()) { // Special condition for 2007/2008 cosmic data
197 0 : dStrip = fgkCosmic2008StripShifts[ddl][ad-1];
198 0 : if (TMath::Abs(dStrip) > 1.5){
199 0 : AliWarning(Form("Indexing error in Cosmic calibration: ddl = %d, dStrip %f\n",ddl,dStrip));
200 : dStrip = 0;
201 0 : }
202 : }
203 :
204 0 : for( int side=0; side<=1; side++ ){
205 :
206 0 : Int_t lab[3]={-2,-2,-2};
207 : Float_t q = 0.;
208 : Float_t y = 0.;
209 : Int_t nDigits = 0;
210 : Int_t ostrip = -2;
211 : Bool_t snFlag = 0;
212 :
213 0 : Int_t n = nStrips[adc][side];
214 0 : for( int istr = 0; istr<n+1; istr++ ){
215 :
216 : bool stripOK = 1;
217 : Int_t strip=0, signal=0;
218 : Float_t noise=1, gain=0;
219 :
220 0 : if( istr<n ){
221 0 : strip = strips[adc][side][istr][0];
222 0 : signal = strips[adc][side][istr][1];
223 :
224 : //cout<<"strip "<<adc<<" / "<<side<<": "<<strip<<endl;
225 :
226 0 : if( cal ){
227 0 : noise = side ?cal->GetNoiseN(strip) :cal->GetNoiseP(strip);
228 0 : gain = side ?cal->GetGainN(strip) :cal->GetGainP(strip);
229 0 : stripOK = ( noise>=1. && signal>=3*noise
230 : //&& !cal->IsPChannelBad(strip)
231 : );
232 0 : }
233 : } else stripOK = 0; // end of data
234 :
235 0 : bool newCluster = ( TMath::Abs(strip-ostrip)>1 || !stripOK );
236 :
237 0 : if( newCluster ){
238 :
239 : //* Store the previous cluster
240 :
241 0 : if( nDigits>0 && q>0 && snFlag ){
242 :
243 0 : if (nClusters1D[side] >= kMaxADCClusters-1 ) {
244 0 : AliWarning("HLT ClustersFinderSSD: Too many 1D clusters !");
245 0 : }else {
246 :
247 0 : Ali1Dcluster &cluster = clusters1D[side][nClusters1D[side]++];
248 0 : cluster.SetY( y / q + dStrip);
249 0 : cluster.SetQ(q);
250 0 : cluster.SetNd(nDigits);
251 0 : cluster.SetLabels(lab);
252 : //cout<<"cluster 1D side "<<side<<": y= "<<y<<" q= "<<q<<" d="<<dStrip<<" Y="<<cluster.GetY()<<endl;
253 : //Split suspiciously big cluster
254 :
255 0 : if( fRecoParam->GetUseUnfoldingInClusterFinderSSD()
256 0 : && nDigits > 4 && nDigits < 25
257 : ){
258 0 : cluster.SetY(y/q + dStrip - 0.25*nDigits);
259 0 : cluster.SetQ(0.5*q);
260 0 : Ali1Dcluster& cluster2 = clusters1D[side][nClusters1D[side]++];
261 0 : cluster2.SetY(y/q + dStrip + 0.25*nDigits);
262 0 : cluster2.SetQ(0.5*q);
263 0 : cluster2.SetNd(nDigits);
264 0 : cluster2.SetLabels(lab);
265 0 : } // unfolding is on
266 : }
267 : }
268 : y = q = 0.;
269 : nDigits = 0;
270 : snFlag = 0;
271 :
272 0 : } //* End store the previous cluster
273 :
274 0 : if( stripOK ){ // add new signal to the cluster
275 0 : signal = (Int_t) ( signal * gain ); // signal is corrected for gain
276 0 : if( signal>fgkThreshold*noise) snFlag = 1;
277 0 : if( cal ) signal = (Int_t) cal->ADCToKeV( signal ); // signal is converted in KeV
278 0 : q += signal; // add digit to current cluster
279 0 : y += strip * signal;
280 0 : nDigits++;
281 : //nstat[side]++;
282 : ostrip = strip;
283 : //cout<<"strip "<<adc<<" / "<<side<<": "<<strip<<" / "<<signal<<" stored"<<endl;
284 :
285 0 : }
286 : } //* end loop over strips
287 :
288 0 : } //* end loop over ADC sides
289 :
290 :
291 : //* 2D clusterfinder
292 :
293 0 : if( nClusters1D[0] && nClusters1D[1] && fModule>=0 ){
294 0 : FindClustersSSD( clusters1D[0], nClusters1D[0], clusters1D[1], nClusters1D[1], clusters);
295 0 : }
296 : //cout<<"SG: "<<ddl<<" "<<ad<<" "<<adc<<": strips "<<nstat[0]<<"+"<<nstat[1]<<", clusters 1D= "<<nClusters1D[0]<<" + "<<nClusters1D[1]<<", 2D= "<<clusters.size()<<endl;
297 :
298 0 : }//* end loop over adc
299 :
300 0 : }//* end of reconstruction of previous module
301 :
302 0 : if( newModule ){
303 :
304 : //*
305 : //* Clean up arrays and set new module
306 : //*
307 :
308 0 : for( int i=0; i<kNADC; i++ ){
309 0 : nStrips[i][0] = 0;
310 0 : nStrips[i][1] = 0;
311 : }
312 : ddl = newDDL;
313 : ad = newAD;
314 0 : }
315 :
316 :
317 : //*
318 : //* Exit main loop when there is no more input
319 : //*
320 :
321 0 : if( !next ) break;
322 :
323 : //*
324 : //* Fill the current strip information
325 : //*
326 :
327 0 : Int_t adc = fRawStream->GetADC();
328 0 : if( adc<0 || adc>=kNADC+2 || (adc>5&&adc<8) ){
329 0 : AliWarning(Form("HLT ClustersFinderSSD: Corrupted data: wrong adc number (%d)", adc));
330 0 : continue;
331 : }
332 :
333 0 : if( adc>7 ) adc-= 2; // shift ADC numbers 8-13 to 6-11
334 :
335 0 : Bool_t side = fRawStream->GetSideFlag();
336 0 : Int_t strip = fRawStream->GetStrip();
337 0 : Int_t signal = fRawStream->GetSignal();
338 :
339 : //cout<<"SSD: "<<ddl<<" "<<ad<<" "<<adc<<" "<<side<<" "<<strip<<" : "<<signal<<endl;
340 :
341 0 : if( strip<0 || strip>767 ){
342 0 : continue;
343 : }
344 :
345 0 : int &n = nStrips[adc][side];
346 0 : if( n >0 ){
347 0 : Int_t oldStrip = strips[adc][side][n-1][0];
348 0 : if( strip==oldStrip ){
349 0 : AliWarning(Form("HLT ClustersFinderSSD: Corrupted data: duplicated signal: ddl %d ad %d adc %d, side %d, strip %d",
350 : ddl, ad, adc, side, strip ));
351 0 : continue;
352 : }
353 0 : }
354 0 : strips[adc][side][n][0] = strip;
355 0 : strips[adc][side][n][1] = signal;
356 0 : n++;
357 :
358 : //cout<<"SSD: "<<fRawStream->GetDDL()<<" "<<fRawStream->GetAD()<<" "
359 : //<<fRawStream->GetADC()<<" "<<fRawStream->GetSideFlag()<<" "<<((int)fRawStream->GetStrip())<<" "<<strip<<" : "<<fRawStream->GetSignal()<<endl;
360 :
361 0 : } //* End main loop over the input
362 :
363 0 : }
364 :
365 :
366 : void AliHLTITSClusterFinderSSD::
367 : FindClustersSSD(Ali1Dcluster* neg, Int_t nn,
368 : Ali1Dcluster* pos, Int_t np,
369 : std::vector<AliITSRecPoint> &clusters) {
370 : //------------------------------------------------------------
371 : // Actual SSD cluster finder
372 : //------------------------------------------------------------
373 :
374 0 : const TGeoHMatrix *mT2L=AliITSgeomTGeo::GetTracking2LocalMatrix(fModule);
375 :
376 0 : AliITSsegmentationSSD *seg = dynamic_cast<AliITSsegmentationSSD*>(fDetTypeRec->GetSegmentationModel(2));
377 0 : if(!mT2L || !seg){
378 0 : AliError(Form("HLT ClustersFinderSSD: null pointer from TGeo"));
379 0 : return;
380 : }
381 0 : if (fModule>fLastSSD1)
382 0 : seg->SetLayer(6);
383 : else
384 0 : seg->SetLayer(5);
385 :
386 0 : Float_t hwSSD = seg->Dx()*1e-4/2;
387 0 : Float_t hlSSD = seg->Dz()*1e-4/2;
388 :
389 0 : Int_t idet=fNdet[fModule];
390 : Int_t ncl=0;
391 :
392 : //
393 0 : Int_t *cnegative = new Int_t[np];
394 0 : Int_t *cused1 = new Int_t[np];
395 0 : Int_t *negativepair = new Int_t[10*np];
396 0 : Int_t *cpositive = new Int_t[nn];
397 0 : Int_t *cused2 = new Int_t[nn];
398 0 : Int_t *positivepair = new Int_t[10*nn];
399 0 : for (Int_t i=0;i<np;i++) {cnegative[i]=0; cused1[i]=0;}
400 0 : for (Int_t i=0;i<nn;i++) {cpositive[i]=0; cused2[i]=0;}
401 0 : for (Int_t i=0;i<10*np;i++) {negativepair[i]=0;}
402 0 : for (Int_t i=0;i<10*nn;i++) {positivepair[i]=0;}
403 :
404 0 : if ((np*nn) > fgPairsSize) {
405 :
406 0 : if (fgPairs) delete [] fgPairs;
407 0 : fgPairsSize = 4*np*nn;
408 0 : fgPairs = new Short_t[fgPairsSize];
409 0 : }
410 0 : memset(fgPairs,0,sizeof(Short_t)*np*nn);
411 :
412 : //
413 : // find available pairs
414 : //
415 0 : for (Int_t i=0; i<np; i++) {
416 0 : Float_t yp=pos[i].GetY();
417 0 : if ( (pos[i].GetQ()>0) && (pos[i].GetQ()<3) ) continue;
418 0 : for (Int_t j=0; j<nn; j++) {
419 0 : if ( (neg[j].GetQ()>0) && (neg[j].GetQ()<3) ) continue;
420 0 : Float_t yn=neg[j].GetY();
421 :
422 0 : Float_t xt, zt;
423 0 : seg->GetPadCxz(yn, yp, xt, zt);
424 : //cout<<yn<<" "<<yp<<" "<<xt<<" "<<zt<<endl;
425 :
426 0 : if (TMath::Abs(xt)<hwSSD)
427 0 : if (TMath::Abs(zt)<hlSSD) {
428 0 : Int_t in = i*10+cnegative[i];
429 0 : Int_t ip = j*10+cpositive[j];
430 0 : if ((in < 10*np) && (ip < 10*nn)) {
431 0 : negativepair[in] =j; //index
432 0 : positivepair[ip] =i;
433 0 : cnegative[i]++; //counters
434 0 : cpositive[j]++;
435 0 : fgPairs[i*nn+j]=100;
436 0 : }
437 : else
438 0 : AliError(Form("Index out of range: ip=%d, in=%d",ip,in));
439 0 : }
440 0 : }
441 0 : }
442 :
443 :
444 : //
445 0 : Float_t lp[6];
446 0 : Int_t milab[10];
447 : Double_t ratio;
448 :
449 :
450 0 : if(fRecoParam->GetUseChargeMatchingInClusterFinderSSD()==kTRUE) {
451 :
452 :
453 : //
454 : // sign gold tracks
455 : //
456 0 : for (Int_t ip=0;ip<np;ip++){
457 :
458 : Float_t xbest=1000,zbest=1000,qbest=0;
459 : //
460 : // select gold clusters
461 0 : if ( (cnegative[ip]==1) && cpositive[negativepair[10*ip]]==1){
462 :
463 0 : Float_t yp=pos[ip].GetY();
464 0 : Int_t j = negativepair[10*ip];
465 :
466 0 : if( (pos[ip].GetQ()==0) && (neg[j].GetQ() ==0) ) {
467 : // both bad, hence continue;
468 : // mark both as used (to avoid recover at the end)
469 0 : cused1[ip]++;
470 0 : cused2[j]++;
471 0 : continue;
472 : }
473 :
474 0 : ratio = (pos[ip].GetQ()-neg[j].GetQ())/(pos[ip].GetQ()+neg[j].GetQ());
475 : //cout<<"ratio="<<ratio<<endl;
476 :
477 : // charge matching (note that if posQ or negQ is 0 -> ratio=1 and the following condition is met
478 0 : if (TMath::Abs(ratio)>0.2) continue; // note: 0.2=3xsigma_ratio calculated in cosmics tests
479 :
480 : //
481 0 : Float_t yn=neg[j].GetY();
482 :
483 0 : Float_t xt, zt;
484 0 : seg->GetPadCxz(yn, yp, xt, zt);
485 :
486 0 : xbest=xt; zbest=zt;
487 :
488 :
489 0 : qbest=0.5*(pos[ip].GetQ()+neg[j].GetQ());
490 0 : 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
491 :
492 : {
493 0 : Double_t loc[3]={xbest,0.,zbest},trk[3]={0.,0.,0.};
494 0 : mT2L->MasterToLocal(loc,trk);
495 0 : lp[0]=trk[1];
496 0 : lp[1]=trk[2];
497 0 : }
498 0 : lp[4]=qbest; //Q
499 0 : for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2;
500 0 : for (Int_t ilab=0;ilab<3;ilab++){
501 0 : milab[ilab] = pos[ip].GetLabel(ilab);
502 0 : milab[ilab+3] = neg[j].GetLabel(ilab);
503 : }
504 : //
505 0 : CheckLabels2(milab);
506 0 : milab[3]=(((ip<<10) + j)<<10) + idet; // pos|neg|det
507 0 : Int_t info[3] = {pos[ip].GetNd(),neg[j].GetNd(),fNlayer[fModule]};
508 :
509 0 : lp[2]=0.0022*0.0022; //SigmaY2
510 0 : lp[3]=0.110*0.110; //SigmaZ2
511 : // out-of-diagonal element of covariance matrix
512 0 : if( (info[0]==1) && (info[1]==1) ) lp[5]=-0.00012;
513 0 : else if ( (info[0]>1) && (info[1]>1) ) {
514 0 : lp[2]=0.0016*0.0016; //SigmaY2
515 0 : lp[3]=0.08*0.08; //SigmaZ2
516 0 : lp[5]=-0.00006;
517 0 : }
518 : else {
519 0 : lp[3]=0.093*0.093;
520 0 : if (info[0]==1) { lp[5]=-0.00014;}
521 0 : else { lp[2]=0.0017*0.0017; lp[5]=-0.00004;}
522 : }
523 0 : AliITSRecPoint cl2(milab,lp,info);
524 0 : cl2.SetChargeRatio(ratio);
525 0 : cl2.SetType(1);
526 :
527 0 : fgPairs[ip*nn+j]=1;
528 0 : if ((pos[ip].GetNd()+neg[j].GetNd())>6){ //multi cluster
529 0 : cl2.SetType(2);
530 0 : fgPairs[ip*nn+j]=2;
531 0 : }
532 :
533 0 : if(pos[ip].GetQ()==0) cl2.SetType(3);
534 0 : if(neg[j].GetQ()==0) cl2.SetType(4);
535 :
536 0 : cused1[ip]++;
537 0 : cused2[j]++;
538 :
539 0 : clusters.push_back(cl2);
540 :
541 0 : ncl++;
542 0 : }
543 0 : }
544 :
545 0 : for (Int_t ip=0;ip<np;ip++){
546 : Float_t xbest=1000,zbest=1000,qbest=0;
547 : //
548 : //
549 : // select "silber" cluster
550 0 : if ( cnegative[ip]==1 && cpositive[negativepair[10*ip]]==2){
551 : Int_t in = negativepair[10*ip];
552 0 : Int_t ip2 = positivepair[10*in];
553 0 : if (ip2==ip) ip2 = positivepair[10*in+1];
554 0 : Float_t pcharge = pos[ip].GetQ()+pos[ip2].GetQ();
555 :
556 :
557 :
558 0 : ratio = (pcharge-neg[in].GetQ())/(pcharge+neg[in].GetQ());
559 0 : if ( (TMath::Abs(ratio)<0.2) && (pcharge!=0) ) {
560 : //if ( (TMath::Abs(pcharge-neg[in].GetQ())<30) && (pcharge!=0) ) { //
561 :
562 : //
563 : // add first pair
564 0 : if ( (fgPairs[ip*nn+in]==100)&&(pos[ip].GetQ() ) ) { //
565 :
566 0 : Float_t yp=pos[ip].GetY();
567 0 : Float_t yn=neg[in].GetY();
568 :
569 0 : Float_t xt, zt;
570 0 : seg->GetPadCxz(yn, yp, xt, zt);
571 :
572 0 : xbest=xt; zbest=zt;
573 :
574 0 : qbest =pos[ip].GetQ();
575 0 : Double_t loc[3]={xbest,0.,zbest},trk[3]={0.,0.,0.};
576 0 : mT2L->MasterToLocal(loc,trk);
577 0 : lp[0]=trk[1];
578 0 : lp[1]=trk[2];
579 :
580 0 : lp[4]=qbest; //Q
581 0 : for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2;
582 0 : for (Int_t ilab=0;ilab<3;ilab++){
583 0 : milab[ilab] = pos[ip].GetLabel(ilab);
584 0 : milab[ilab+3] = neg[in].GetLabel(ilab);
585 : }
586 : //
587 0 : CheckLabels2(milab);
588 0 : ratio = (pos[ip].GetQ()-neg[in].GetQ())/(pos[ip].GetQ()+neg[in].GetQ());
589 0 : milab[3]=(((ip<<10) + in)<<10) + idet; // pos|neg|det
590 0 : Int_t info[3] = {pos[ip].GetNd(),neg[in].GetNd(),fNlayer[fModule]};
591 :
592 0 : lp[2]=0.0022*0.0022; //SigmaY2
593 0 : lp[3]=0.110*0.110; //SigmaZ2
594 : // out-of-diagonal element of covariance matrix
595 0 : if( (info[0]==1) && (info[1]==1) ) lp[5]=-0.00012;
596 0 : else if ( (info[0]>1) && (info[1]>1) ) {
597 0 : lp[2]=0.0016*0.0016; //SigmaY2
598 0 : lp[3]=0.08*0.08; //SigmaZ2
599 0 : lp[5]=-0.00006;
600 0 : }
601 : else {
602 0 : lp[3]=0.093*0.093;
603 0 : if (info[0]==1) { lp[5]=-0.00014;}
604 0 : else { lp[2]=0.0017*0.0017; lp[5]=-0.00004;}
605 : }
606 :
607 0 : AliITSRecPoint cl2(milab,lp,info);
608 0 : cl2.SetChargeRatio(ratio);
609 0 : cl2.SetType(5);
610 0 : fgPairs[ip*nn+in] = 5;
611 0 : if ((pos[ip].GetNd()+neg[in].GetNd())>6){ //multi cluster
612 0 : cl2.SetType(6);
613 0 : fgPairs[ip*nn+in] = 6;
614 0 : }
615 0 : clusters.push_back(cl2);
616 0 : ncl++;
617 0 : }
618 :
619 :
620 : //
621 : // add second pair
622 :
623 : // if (!(cused1[ip2] || cused2[in])){ //
624 0 : if ( (fgPairs[ip2*nn+in]==100) && (pos[ip2].GetQ()) ) {
625 :
626 0 : Float_t yp=pos[ip2].GetY();
627 0 : Float_t yn=neg[in].GetY();
628 :
629 0 : Float_t xt, zt;
630 0 : seg->GetPadCxz(yn, yp, xt, zt);
631 :
632 0 : xbest=xt; zbest=zt;
633 :
634 0 : qbest =pos[ip2].GetQ();
635 :
636 0 : Double_t loc[3]={xbest,0.,zbest},trk[3]={0.,0.,0.};
637 0 : mT2L->MasterToLocal(loc,trk);
638 0 : lp[0]=trk[1];
639 0 : lp[1]=trk[2];
640 :
641 0 : lp[4]=qbest; //Q
642 0 : for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2;
643 0 : for (Int_t ilab=0;ilab<3;ilab++){
644 0 : milab[ilab] = pos[ip2].GetLabel(ilab);
645 0 : milab[ilab+3] = neg[in].GetLabel(ilab);
646 : }
647 : //
648 0 : CheckLabels2(milab);
649 0 : ratio = (pos[ip2].GetQ()-neg[in].GetQ())/(pos[ip2].GetQ()+neg[in].GetQ());
650 0 : milab[3]=(((ip2<<10) + in)<<10) + idet; // pos|neg|det
651 0 : Int_t info[3] = {pos[ip2].GetNd(),neg[in].GetNd(),fNlayer[fModule]};
652 :
653 0 : lp[2]=0.0022*0.0022; //SigmaY2
654 0 : lp[3]=0.110*0.110; //SigmaZ2
655 : // out-of-diagonal element of covariance matrix
656 0 : if( (info[0]==1) && (info[1]==1) ) lp[5]=-0.00012;
657 0 : else if ( (info[0]>1) && (info[1]>1) ) {
658 0 : lp[2]=0.0016*0.0016; //SigmaY2
659 0 : lp[3]=0.08*0.08; //SigmaZ2
660 0 : lp[5]=-0.00006;
661 0 : }
662 : else {
663 0 : lp[3]=0.093*0.093;
664 0 : if (info[0]==1) { lp[5]=-0.00014;}
665 0 : else { lp[2]=0.0017*0.0017; lp[5]=-0.00004;}
666 : }
667 :
668 0 : AliITSRecPoint cl2(milab,lp,info);
669 0 : cl2.SetChargeRatio(ratio);
670 0 : cl2.SetType(5);
671 0 : fgPairs[ip2*nn+in] =5;
672 0 : if ((pos[ip2].GetNd()+neg[in].GetNd())>6){ //multi cluster
673 0 : cl2.SetType(6);
674 0 : fgPairs[ip2*nn+in] =6;
675 0 : }
676 0 : clusters.push_back(cl2);
677 0 : ncl++;
678 0 : }
679 :
680 0 : cused1[ip]++;
681 0 : cused1[ip2]++;
682 0 : cused2[in]++;
683 :
684 0 : } // charge matching condition
685 :
686 0 : } // 2 Pside cross 1 Nside
687 : } // loop over Pside clusters
688 :
689 :
690 : //
691 0 : for (Int_t jn=0;jn<nn;jn++){
692 0 : if (cused2[jn]) continue;
693 : Float_t xbest=1000,zbest=1000,qbest=0;
694 : // select "silber" cluster
695 0 : if ( cpositive[jn]==1 && cnegative[positivepair[10*jn]]==2){
696 : Int_t ip = positivepair[10*jn];
697 0 : Int_t jn2 = negativepair[10*ip];
698 0 : if (jn2==jn) jn2 = negativepair[10*ip+1];
699 0 : Float_t pcharge = neg[jn].GetQ()+neg[jn2].GetQ();
700 : //
701 :
702 :
703 0 : ratio = (pcharge-pos[ip].GetQ())/(pcharge+pos[ip].GetQ());
704 0 : if ( (TMath::Abs(ratio)<0.2) && (pcharge!=0) ) {
705 :
706 : /*
707 : if ( (TMath::Abs(pcharge-pos[ip].GetQ())<30) && // charge matching
708 : (pcharge!=0) ) { // reject combinations of bad strips
709 : */
710 :
711 :
712 : //
713 : // add first pair
714 : // if (!(cused1[ip]||cused2[jn])){
715 0 : if ( (fgPairs[ip*nn+jn]==100) && (neg[jn].GetQ()) ) { //
716 :
717 0 : Float_t yn=neg[jn].GetY();
718 0 : Float_t yp=pos[ip].GetY();
719 :
720 0 : Float_t xt, zt;
721 0 : seg->GetPadCxz(yn, yp, xt, zt);
722 :
723 0 : xbest=xt; zbest=zt;
724 :
725 0 : qbest =neg[jn].GetQ();
726 :
727 : {
728 0 : Double_t loc[3]={xbest,0.,zbest},trk[3]={0.,0.,0.};
729 0 : mT2L->MasterToLocal(loc,trk);
730 0 : lp[0]=trk[1];
731 0 : lp[1]=trk[2];
732 0 : }
733 :
734 0 : lp[4]=qbest; //Q
735 0 : for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2;
736 0 : for (Int_t ilab=0;ilab<3;ilab++){
737 0 : milab[ilab] = pos[ip].GetLabel(ilab);
738 0 : milab[ilab+3] = neg[jn].GetLabel(ilab);
739 : }
740 : //
741 0 : CheckLabels2(milab);
742 0 : ratio = (pos[ip].GetQ()-neg[jn].GetQ())/(pos[ip].GetQ()+neg[jn].GetQ());
743 0 : milab[3]=(((ip<<10) + jn)<<10) + idet; // pos|neg|det
744 0 : Int_t info[3] = {pos[ip].GetNd(),neg[jn].GetNd(),fNlayer[fModule]};
745 :
746 0 : lp[2]=0.0022*0.0022; //SigmaY2
747 0 : lp[3]=0.110*0.110; //SigmaZ2
748 : // out-of-diagonal element of covariance matrix
749 0 : if( (info[0]==1) && (info[1]==1) ) lp[5]=-0.00012;
750 0 : else if ( (info[0]>1) && (info[1]>1) ) {
751 0 : lp[2]=0.0016*0.0016; //SigmaY2
752 0 : lp[3]=0.08*0.08; //SigmaZ2
753 0 : lp[5]=-0.00006;
754 0 : }
755 : else {
756 0 : lp[3]=0.093*0.093;
757 0 : if (info[0]==1) { lp[5]=-0.00014;}
758 0 : else { lp[2]=0.0017*0.0017; lp[5]=-0.00004;}
759 : }
760 :
761 0 : AliITSRecPoint cl2(milab,lp,info);
762 0 : cl2.SetChargeRatio(ratio);
763 0 : cl2.SetType(7);
764 0 : fgPairs[ip*nn+jn] =7;
765 0 : if ((pos[ip].GetNd()+neg[jn].GetNd())>6){ //multi cluster
766 0 : cl2.SetType(8);
767 0 : fgPairs[ip*nn+jn]=8;
768 0 : }
769 0 : clusters.push_back(cl2);
770 :
771 0 : ncl++;
772 0 : }
773 : //
774 : // add second pair
775 : // if (!(cused1[ip]||cused2[jn2])){
776 0 : if ( (fgPairs[ip*nn+jn2]==100)&&(neg[jn2].GetQ() ) ) { //
777 :
778 0 : Float_t yn=neg[jn2].GetY();
779 0 : Double_t yp=pos[ip].GetY();
780 :
781 0 : Float_t xt, zt;
782 0 : seg->GetPadCxz(yn, yp, xt, zt);
783 :
784 0 : xbest=xt; zbest=zt;
785 :
786 0 : qbest =neg[jn2].GetQ();
787 :
788 : {
789 0 : Double_t loc[3]={xbest,0.,zbest},trk[3]={0.,0.,0.};
790 0 : mT2L->MasterToLocal(loc,trk);
791 0 : lp[0]=trk[1];
792 0 : lp[1]=trk[2];
793 0 : }
794 :
795 0 : lp[4]=qbest; //Q
796 0 : for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2;
797 0 : for (Int_t ilab=0;ilab<3;ilab++){
798 0 : milab[ilab] = pos[ip].GetLabel(ilab);
799 0 : milab[ilab+3] = neg[jn2].GetLabel(ilab);
800 : }
801 : //
802 0 : CheckLabels2(milab);
803 0 : ratio = (pos[ip].GetQ()-neg[jn2].GetQ())/(pos[ip].GetQ()+neg[jn2].GetQ());
804 0 : milab[3]=(((ip<<10) + jn2)<<10) + idet; // pos|neg|det
805 0 : Int_t info[3] = {pos[ip].GetNd(),neg[jn2].GetNd(),fNlayer[fModule]};
806 :
807 0 : lp[2]=0.0022*0.0022; //SigmaY2
808 0 : lp[3]=0.110*0.110; //SigmaZ2
809 : // out-of-diagonal element of covariance matrix
810 0 : if( (info[0]==1) && (info[1]==1) ) lp[5]=-0.00012;
811 0 : else if ( (info[0]>1) && (info[1]>1) ) {
812 0 : lp[2]=0.0016*0.0016; //SigmaY2
813 0 : lp[3]=0.08*0.08; //SigmaZ2
814 0 : lp[5]=-0.00006;
815 0 : }
816 : else {
817 0 : lp[3]=0.093*0.093;
818 0 : if (info[0]==1) { lp[5]=-0.00014;}
819 0 : else { lp[2]=0.0017*0.0017; lp[5]=-0.00004;}
820 : }
821 :
822 0 : AliITSRecPoint cl2(milab,lp,info);
823 0 : cl2.SetChargeRatio(ratio);
824 0 : fgPairs[ip*nn+jn2]=7;
825 0 : cl2.SetType(7);
826 0 : if ((pos[ip].GetNd()+neg[jn2].GetNd())>6){ //multi cluster
827 0 : cl2.SetType(8);
828 0 : fgPairs[ip*nn+jn2]=8;
829 0 : }
830 0 : clusters.push_back( cl2 );
831 0 : ncl++;
832 0 : }
833 0 : cused1[ip]++;
834 0 : cused2[jn]++;
835 0 : cused2[jn2]++;
836 :
837 0 : } // charge matching condition
838 :
839 0 : } // 2 Nside cross 1 Pside
840 0 : } // loop over Pside clusters
841 :
842 :
843 :
844 0 : for (Int_t ip=0;ip<np;ip++){
845 :
846 0 : if(cused1[ip]) continue;
847 :
848 :
849 : Float_t xbest=1000,zbest=1000,qbest=0;
850 : //
851 : // 2x2 clusters
852 : //
853 0 : if ( (cnegative[ip]==2) && cpositive[negativepair[10*ip]]==2){
854 : Float_t minchargediff =4.;
855 : //Float_t minchargeratio =0.2;
856 :
857 : Int_t j=-1;
858 0 : for (Int_t di=0;di<cnegative[ip];di++){
859 0 : Int_t jc = negativepair[ip*10+di];
860 0 : Float_t chargedif = pos[ip].GetQ()-neg[jc].GetQ();
861 0 : ratio = (pos[ip].GetQ()-neg[jc].GetQ())/(pos[ip].GetQ()+neg[jc].GetQ());
862 : //if (TMath::Abs(chargedif)<minchargediff){
863 0 : if (TMath::Abs(ratio)<0.2){
864 : j =jc;
865 0 : minchargediff = TMath::Abs(chargedif);
866 : //minchargeratio = TMath::Abs(ratio);
867 0 : }
868 : }
869 0 : if (j<0) continue; // not proper cluster
870 :
871 :
872 : Int_t count =0;
873 0 : for (Int_t di=0;di<cnegative[ip];di++){
874 0 : Int_t jc = negativepair[ip*10+di];
875 0 : Float_t chargedif = pos[ip].GetQ()-neg[jc].GetQ();
876 0 : if (TMath::Abs(chargedif)<minchargediff+3.) count++;
877 : }
878 0 : if (count>1) continue; // more than one "proper" cluster for positive
879 : //
880 :
881 : count =0;
882 0 : for (Int_t dj=0;dj<cpositive[j];dj++){
883 0 : Int_t ic = positivepair[j*10+dj];
884 0 : Float_t chargedif = pos[ic].GetQ()-neg[j].GetQ();
885 0 : if (TMath::Abs(chargedif)<minchargediff+3.) count++;
886 : }
887 0 : if (count>1) continue; // more than one "proper" cluster for negative
888 :
889 : Int_t jp = 0;
890 :
891 : count =0;
892 0 : for (Int_t dj=0;dj<cnegative[jp];dj++){
893 0 : Int_t ic = positivepair[jp*10+dj];
894 0 : Float_t chargedif = pos[ic].GetQ()-neg[jp].GetQ();
895 0 : if (TMath::Abs(chargedif)<minchargediff+4.) count++;
896 : }
897 0 : if (count>1) continue;
898 0 : if (fgPairs[ip*nn+j]<100) continue;
899 : //
900 :
901 :
902 :
903 : //almost gold clusters
904 0 : Float_t yp=pos[ip].GetY();
905 0 : Float_t yn=neg[j].GetY();
906 0 : Float_t xt, zt;
907 0 : seg->GetPadCxz(yn, yp, xt, zt);
908 0 : xbest=xt; zbest=zt;
909 0 : qbest=0.5*(pos[ip].GetQ()+neg[j].GetQ());
910 : {
911 0 : Double_t loc[3]={xbest,0.,zbest},trk[3]={0.,0.,0.};
912 0 : mT2L->MasterToLocal(loc,trk);
913 0 : lp[0]=trk[1];
914 0 : lp[1]=trk[2];
915 0 : }
916 0 : lp[4]=qbest; //Q
917 0 : for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2;
918 0 : for (Int_t ilab=0;ilab<3;ilab++){
919 0 : milab[ilab] = pos[ip].GetLabel(ilab);
920 0 : milab[ilab+3] = neg[j].GetLabel(ilab);
921 : }
922 : //
923 0 : CheckLabels2(milab);
924 0 : if ((neg[j].GetQ()==0)&&(pos[ip].GetQ()==0)) continue; // reject crosses of bad strips!!
925 0 : ratio = (pos[ip].GetQ()-neg[j].GetQ())/(pos[ip].GetQ()+neg[j].GetQ());
926 0 : milab[3]=(((ip<<10) + j)<<10) + idet; // pos|neg|det
927 0 : Int_t info[3] = {pos[ip].GetNd(),neg[j].GetNd(),fNlayer[fModule]};
928 :
929 0 : lp[2]=0.0022*0.0022; //SigmaY2
930 0 : lp[3]=0.110*0.110; //SigmaZ2
931 : // out-of-diagonal element of covariance matrix
932 0 : if( (info[0]==1) && (info[1]==1) ) lp[5]=-0.00012;
933 0 : else if ( (info[0]>1) && (info[1]>1) ) {
934 0 : lp[2]=0.0016*0.0016; //SigmaY2
935 0 : lp[3]=0.08*0.08; //SigmaZ2
936 0 : lp[5]=-0.00006;
937 0 : }
938 : else {
939 0 : lp[3]=0.093*0.093;
940 0 : if (info[0]==1) { lp[5]=-0.00014;}
941 0 : else { lp[2]=0.0017*0.0017; lp[5]=-0.00004;}
942 : }
943 :
944 0 : AliITSRecPoint cl2(milab,lp,info);
945 0 : cl2.SetChargeRatio(ratio);
946 0 : cl2.SetType(10);
947 0 : fgPairs[ip*nn+j]=10;
948 0 : if ((pos[ip].GetNd()+neg[j].GetNd())>6){ //multi cluster
949 0 : cl2.SetType(11);
950 0 : fgPairs[ip*nn+j]=11;
951 0 : }
952 0 : cused1[ip]++;
953 0 : cused2[j]++;
954 :
955 0 : clusters.push_back(cl2);
956 0 : ncl++;
957 :
958 0 : } // 2X2
959 0 : } // loop over Pside 1Dclusters
960 :
961 :
962 0 : for (Int_t ip=0;ip<np;ip++){
963 :
964 0 : if(cused1[ip]) continue;
965 :
966 :
967 : Float_t xbest=1000,zbest=1000,qbest=0;
968 : //
969 : // manyxmany clusters
970 : //
971 0 : if ( (cnegative[ip]<5) && cpositive[negativepair[10*ip]]<5){
972 : Float_t minchargediff =4.;
973 : Int_t j=-1;
974 0 : for (Int_t di=0;di<cnegative[ip];di++){
975 0 : Int_t jc = negativepair[ip*10+di];
976 0 : Float_t chargedif = pos[ip].GetQ()-neg[jc].GetQ();
977 0 : if (TMath::Abs(chargedif)<minchargediff){
978 : j =jc;
979 0 : minchargediff = TMath::Abs(chargedif);
980 0 : }
981 : }
982 0 : if (j<0) continue; // not proper cluster
983 :
984 : Int_t count =0;
985 0 : for (Int_t di=0;di<cnegative[ip];di++){
986 0 : Int_t jc = negativepair[ip*10+di];
987 0 : Float_t chargedif = pos[ip].GetQ()-neg[jc].GetQ();
988 0 : if (TMath::Abs(chargedif)<minchargediff+3.) count++;
989 : }
990 0 : if (count>1) continue; // more than one "proper" cluster for positive
991 : //
992 :
993 : count =0;
994 0 : for (Int_t dj=0;dj<cpositive[j];dj++){
995 0 : Int_t ic = positivepair[j*10+dj];
996 0 : Float_t chargedif = pos[ic].GetQ()-neg[j].GetQ();
997 0 : if (TMath::Abs(chargedif)<minchargediff+3.) count++;
998 : }
999 0 : if (count>1) continue; // more than one "proper" cluster for negative
1000 :
1001 : Int_t jp = 0;
1002 :
1003 : count =0;
1004 0 : for (Int_t dj=0;dj<cnegative[jp];dj++){
1005 0 : Int_t ic = positivepair[jp*10+dj];
1006 0 : Float_t chargedif = pos[ic].GetQ()-neg[jp].GetQ();
1007 0 : if (TMath::Abs(chargedif)<minchargediff+4.) count++;
1008 : }
1009 0 : if (count>1) continue;
1010 0 : if (fgPairs[ip*nn+j]<100) continue;
1011 : //
1012 :
1013 : //almost gold clusters
1014 0 : Float_t yp=pos[ip].GetY();
1015 0 : Float_t yn=neg[j].GetY();
1016 :
1017 :
1018 0 : Float_t xt, zt;
1019 0 : seg->GetPadCxz(yn, yp, xt, zt);
1020 :
1021 0 : xbest=xt; zbest=zt;
1022 :
1023 0 : qbest=0.5*(pos[ip].GetQ()+neg[j].GetQ());
1024 :
1025 : {
1026 0 : Double_t loc[3]={xbest,0.,zbest},trk[3]={0.,0.,0.};
1027 0 : mT2L->MasterToLocal(loc,trk);
1028 0 : lp[0]=trk[1];
1029 0 : lp[1]=trk[2];
1030 0 : }
1031 0 : lp[4]=qbest; //Q
1032 0 : for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2;
1033 0 : for (Int_t ilab=0;ilab<3;ilab++){
1034 0 : milab[ilab] = pos[ip].GetLabel(ilab);
1035 0 : milab[ilab+3] = neg[j].GetLabel(ilab);
1036 : }
1037 : //
1038 0 : CheckLabels2(milab);
1039 0 : if ((neg[j].GetQ()==0)&&(pos[ip].GetQ()==0)) continue; // reject crosses of bad strips!!
1040 0 : ratio = (pos[ip].GetQ()-neg[j].GetQ())/(pos[ip].GetQ()+neg[j].GetQ());
1041 0 : milab[3]=(((ip<<10) + j)<<10) + idet; // pos|neg|det
1042 0 : Int_t info[3] = {pos[ip].GetNd(),neg[j].GetNd(),fNlayer[fModule]};
1043 :
1044 0 : lp[2]=0.0022*0.0022; //SigmaY2
1045 0 : lp[3]=0.110*0.110; //SigmaZ2
1046 : // out-of-diagonal element of covariance matrix
1047 0 : if( (info[0]==1) && (info[1]==1) ) lp[5]=-0.00012;
1048 0 : else if ( (info[0]>1) && (info[1]>1) ) {
1049 0 : lp[2]=0.0016*0.0016; //SigmaY2
1050 0 : lp[3]=0.08*0.08; //SigmaZ2
1051 0 : lp[5]=-0.00006;
1052 0 : }
1053 : else {
1054 0 : lp[3]=0.093*0.093;
1055 0 : if (info[0]==1) { lp[5]=-0.00014;}
1056 0 : else { lp[2]=0.0017*0.0017; lp[5]=-0.00004;}
1057 : }
1058 :
1059 0 : AliITSRecPoint cl2(milab,lp,info);
1060 0 : cl2.SetChargeRatio(ratio);
1061 0 : cl2.SetType(12);
1062 0 : fgPairs[ip*nn+j]=12;
1063 0 : if ((pos[ip].GetNd()+neg[j].GetNd())>6){ //multi cluster
1064 0 : cl2.SetType(13);
1065 0 : fgPairs[ip*nn+j]=13;
1066 0 : }
1067 0 : cused1[ip]++;
1068 0 : cused2[j]++;
1069 0 : clusters.push_back( cl2 );
1070 0 : ncl++;
1071 :
1072 0 : } // manyXmany
1073 0 : } // loop over Pside 1Dclusters
1074 :
1075 0 : } // use charge matching
1076 :
1077 : // recover all the other crosses
1078 : //
1079 0 : for (Int_t i=0; i<np; i++) {
1080 : Float_t xbest=1000,zbest=1000,qbest=0;
1081 0 : Float_t yp=pos[i].GetY();
1082 0 : if ((pos[i].GetQ()>0)&&(pos[i].GetQ()<3)) continue;
1083 0 : for (Int_t j=0; j<nn; j++) {
1084 : // for (Int_t di = 0;di<cpositive[i];di++){
1085 : // Int_t j = negativepair[10*i+di];
1086 0 : if ((neg[j].GetQ()>0)&&(neg[j].GetQ()<3)) continue;
1087 :
1088 0 : if ((neg[j].GetQ()==0)&&(pos[i].GetQ()==0)) continue; // reject crosses of bad strips!!
1089 :
1090 0 : if (cused2[j]||cused1[i]) continue;
1091 0 : if (fgPairs[i*nn+j]>0 &&fgPairs[i*nn+j]<100) continue;
1092 0 : ratio = (pos[i].GetQ()-neg[j].GetQ())/(pos[i].GetQ()+neg[j].GetQ());
1093 0 : Float_t yn=neg[j].GetY();
1094 :
1095 0 : Float_t xt, zt;
1096 0 : seg->GetPadCxz(yn, yp, xt, zt);
1097 :
1098 0 : if (TMath::Abs(xt)<hwSSD)
1099 0 : if (TMath::Abs(zt)<hlSSD) {
1100 0 : xbest=xt; zbest=zt;
1101 :
1102 0 : qbest=0.5*(pos[i].GetQ()+neg[j].GetQ());
1103 :
1104 : {
1105 0 : Double_t loc[3]={xbest,0.,zbest},trk[3]={0.,0.,0.};
1106 0 : mT2L->MasterToLocal(loc,trk);
1107 0 : lp[0]=trk[1];
1108 0 : lp[1]=trk[2];
1109 0 : }
1110 0 : lp[4]=qbest; //Q
1111 0 : for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2;
1112 0 : for (Int_t ilab=0;ilab<3;ilab++){
1113 0 : milab[ilab] = pos[i].GetLabel(ilab);
1114 0 : milab[ilab+3] = neg[j].GetLabel(ilab);
1115 : }
1116 : //
1117 0 : CheckLabels2(milab);
1118 0 : milab[3]=(((i<<10) + j)<<10) + idet; // pos|neg|det
1119 0 : Int_t info[3] = {pos[i].GetNd(),neg[j].GetNd(),fNlayer[fModule]};
1120 :
1121 0 : lp[2]=0.0022*0.0022; //SigmaY2
1122 0 : lp[3]=0.110*0.110; //SigmaZ2
1123 : // out-of-diagonal element of covariance matrix
1124 0 : if( (info[0]==1) && (info[1]==1) ) lp[5]=-0.00012;
1125 0 : else if ( (info[0]>1) && (info[1]>1) ) {
1126 0 : lp[2]=0.0016*0.0016; //SigmaY2
1127 0 : lp[3]=0.08*0.08; //SigmaZ2
1128 0 : lp[5]=-0.00006;
1129 0 : }
1130 : else {
1131 0 : lp[3]=0.093*0.093;
1132 0 : if (info[0]==1) { lp[5]=-0.00014;}
1133 0 : else { lp[2]=0.0017*0.0017; lp[5]=-0.00004;}
1134 : }
1135 :
1136 0 : AliITSRecPoint cl2(milab,lp,info);
1137 0 : cl2.SetChargeRatio(ratio);
1138 0 : cl2.SetType(100+cpositive[j]+cnegative[i]);
1139 :
1140 0 : if(pos[i].GetQ()==0) cl2.SetType(200+cpositive[j]+cnegative[i]);
1141 0 : if(neg[j].GetQ()==0) cl2.SetType(300+cpositive[j]+cnegative[i]);
1142 0 : clusters.push_back( cl2 );
1143 0 : ncl++;
1144 0 : }
1145 0 : }
1146 0 : }
1147 :
1148 :
1149 0 : if(fRecoParam->GetUseBadChannelsInClusterFinderSSD()==kTRUE) {
1150 :
1151 : //---------------------------------------------------------
1152 : // recover crosses of good 1D clusters with bad strips on the other side
1153 : // Note1: at first iteration skip modules with a bad side (or almost), (would produce too many fake!)
1154 : // Note2: for modules with a bad side see below
1155 :
1156 0 : AliITSCalibrationSSD* cal = (AliITSCalibrationSSD*) fDetTypeRec->GetCalibrationModel(fModule);
1157 : Int_t countPbad=0, countNbad=0;
1158 0 : for(Int_t ib=0; ib<768; ib++) {
1159 0 : if(cal->IsPChannelBad(ib)) countPbad++;
1160 0 : if(cal->IsNChannelBad(ib)) countNbad++;
1161 : }
1162 : // AliInfo(Form("module %d has %d P- and %d N-bad strips",fModule,countPbad,countNbad));
1163 :
1164 0 : if( (countPbad<100) && (countNbad<100) ) { // no bad side!!
1165 :
1166 0 : for (Int_t i=0; i<np; i++) { // loop over Nside 1Dclusters with no crosses
1167 0 : if(cnegative[i]) continue; // if intersecting Pside clusters continue;
1168 :
1169 : // for(Int_t ib=0; ib<768; ib++) { // loop over all Pstrips
1170 0 : for(Int_t ib=15; ib<753; ib++) { // loop over all Pstrips
1171 :
1172 0 : if(cal->IsPChannelBad(ib)) { // check if strips is bad
1173 0 : Float_t yN=pos[i].GetY();
1174 0 : Float_t xt, zt;
1175 0 : seg->GetPadCxz(1.*ib, yN, xt, zt);
1176 :
1177 : //----------
1178 : // bad Pstrip is crossing the Nside 1Dcluster -> create recpoint
1179 : //
1180 0 : if ( (TMath::Abs(xt)<hwSSD) && (TMath::Abs(zt)<hlSSD) ) {
1181 0 : Double_t loc[3]={xt,0.,zt},trk[3]={0.,0.,0.};
1182 0 : mT2L->MasterToLocal(loc,trk);
1183 0 : lp[0]=trk[1];
1184 0 : lp[1]=trk[2];
1185 0 : lp[4]=pos[i].GetQ(); //Q
1186 0 : for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2;
1187 0 : for (Int_t ilab=0;ilab<3;ilab++) milab[ilab] = pos[i].GetLabel(ilab);
1188 0 : CheckLabels2(milab);
1189 0 : milab[3]=( (i<<10) << 10 ) + idet; // pos|neg|det
1190 0 : Int_t info[3] = {pos[i].GetNd(),0,fNlayer[fModule]};
1191 :
1192 : // out-of-diagonal element of covariance matrix
1193 0 : if (info[0]==1) lp[5]=0.0065;
1194 0 : else lp[5]=0.0093;
1195 :
1196 0 : lp[2]=0.0022*0.0022; //SigmaY2
1197 0 : lp[3]=0.110*0.110; //SigmaZ2
1198 0 : lp[5]=-0.00012; // out-of-diagonal element of covariance matrix
1199 :
1200 0 : AliITSRecPoint cl2(milab,lp,info);
1201 0 : cl2.SetChargeRatio(1.);
1202 0 : cl2.SetType(50);
1203 0 : clusters.push_back( cl2 );
1204 0 : ncl++;
1205 0 : } // cross is within the detector
1206 : //
1207 : //--------------
1208 :
1209 0 : } // bad Pstrip
1210 :
1211 : } // end loop over Pstrips
1212 :
1213 0 : } // end loop over Nside 1D clusters
1214 :
1215 0 : for (Int_t j=0; j<nn; j++) { // loop over Pside 1D clusters with no crosses
1216 0 : if(cpositive[j]) continue;
1217 :
1218 : // for(Int_t ib=0; ib<768; ib++) { // loop over all Nside strips
1219 0 : for(Int_t ib=15; ib<753; ib++) { // loop over all Nside strips
1220 :
1221 0 : if(cal->IsNChannelBad(ib)) { // check if strip is bad
1222 0 : Float_t yP=neg[j].GetY();
1223 0 : Float_t xt, zt;
1224 0 : seg->GetPadCxz(yP, 1.*ib, xt, zt);
1225 :
1226 : //----------
1227 : // bad Nstrip is crossing the Pside 1Dcluster -> create recpoint
1228 : //
1229 0 : if ( (TMath::Abs(xt)<hwSSD) && (TMath::Abs(zt)<hlSSD) ) {
1230 0 : Double_t loc[3]={xt,0.,zt},trk[3]={0.,0.,0.};
1231 0 : mT2L->MasterToLocal(loc,trk);
1232 0 : lp[0]=trk[1];
1233 0 : lp[1]=trk[2];
1234 0 : lp[4]=neg[j].GetQ(); //Q
1235 0 : for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2;
1236 0 : for (Int_t ilab=0;ilab<3;ilab++) milab[ilab] = neg[j].GetLabel(ilab);
1237 0 : CheckLabels2(milab);
1238 0 : milab[3]=( j << 10 ) + idet; // pos|neg|det
1239 0 : Int_t info[3]={0,(Int_t)neg[j].GetNd(),fNlayer[fModule]};
1240 :
1241 0 : lp[2]=0.0022*0.0022; //SigmaY2
1242 0 : lp[3]=0.110*0.110; //SigmaZ2
1243 0 : lp[5]=-0.00012; // out-of-diagonal element of covariance matrix
1244 :
1245 0 : AliITSRecPoint cl2(milab,lp,info);
1246 0 : cl2.SetChargeRatio(1.);
1247 0 : cl2.SetType(60);
1248 0 : clusters.push_back( cl2 );
1249 0 : ncl++;
1250 0 : } // cross is within the detector
1251 : //
1252 : //--------------
1253 :
1254 0 : } // bad Nstrip
1255 : } // end loop over Nstrips
1256 0 : } // end loop over Pside 1D clusters
1257 :
1258 0 : } // no bad sides
1259 :
1260 : //---------------------------------------------------------
1261 :
1262 0 : else if( (countPbad>700) && (countNbad<100) ) { // bad Pside!!
1263 :
1264 0 : for (Int_t i=0; i<np; i++) { // loop over Nside 1Dclusters with no crosses
1265 0 : if(cnegative[i]) continue; // if intersecting Pside clusters continue;
1266 :
1267 0 : Float_t xt, zt;
1268 0 : Float_t yN=pos[i].GetY();
1269 : Float_t yP=0.;
1270 0 : if (seg->GetLayer()==5) yP = yN + (7.6/1.9);
1271 0 : else yP = yN - (7.6/1.9);
1272 0 : seg->GetPadCxz(yP, yN, xt, zt);
1273 :
1274 0 : if ( (TMath::Abs(xt)<hwSSD) && (TMath::Abs(zt)<hlSSD) ) {
1275 0 : Double_t loc[3]={xt,0.,zt},trk[3]={0.,0.,0.};
1276 0 : mT2L->MasterToLocal(loc,trk);
1277 0 : lp[0]=trk[1];
1278 0 : lp[1]=trk[2];
1279 0 : lp[4]=pos[i].GetQ(); //Q
1280 0 : for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2;
1281 0 : for (Int_t ilab=0;ilab<3;ilab++) milab[ilab] = pos[i].GetLabel(ilab);
1282 0 : CheckLabels2(milab);
1283 0 : milab[3]=( (i<<10) << 10 ) + idet; // pos|neg|det
1284 0 : Int_t info[3] = {(Int_t)pos[i].GetNd(),0,fNlayer[fModule]};
1285 :
1286 0 : lp[2]=0.031*0.031; //SigmaY2
1287 0 : lp[3]=1.15*1.15; //SigmaZ2
1288 0 : lp[5]=-0.036;
1289 :
1290 0 : AliITSRecPoint cl2(milab,lp,info);
1291 0 : cl2.SetChargeRatio(1.);
1292 0 : cl2.SetType(70);
1293 0 : clusters.push_back( cl2 );
1294 0 : ncl++;
1295 0 : } // cross is within the detector
1296 : //
1297 : //--------------
1298 :
1299 0 : } // end loop over Nside 1D clusters
1300 :
1301 0 : } // bad Pside module
1302 :
1303 0 : else if( (countNbad>700) && (countPbad<100) ) { // bad Nside!!
1304 :
1305 0 : for (Int_t j=0; j<nn; j++) { // loop over Pside 1D clusters with no crosses
1306 0 : if(cpositive[j]) continue;
1307 :
1308 0 : Float_t xt, zt;
1309 0 : Float_t yP=neg[j].GetY();
1310 : Float_t yN=0.;
1311 0 : if (seg->GetLayer()==5) yN = yP - (7.6/1.9);
1312 0 : else yN = yP + (7.6/1.9);
1313 0 : seg->GetPadCxz(yP, yN, xt, zt);
1314 :
1315 0 : if ( (TMath::Abs(xt)<hwSSD) && (TMath::Abs(zt)<hlSSD) ) {
1316 0 : Double_t loc[3]={xt,0.,zt},trk[3]={0.,0.,0.};
1317 0 : mT2L->MasterToLocal(loc,trk);
1318 0 : lp[0]=trk[1];
1319 0 : lp[1]=trk[2];
1320 0 : lp[4]=neg[j].GetQ(); //Q
1321 0 : for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2;
1322 0 : for (Int_t ilab=0;ilab<3;ilab++) milab[ilab] = neg[j].GetLabel(ilab);
1323 0 : CheckLabels2(milab);
1324 0 : milab[3]=( j << 10 ) + idet; // pos|neg|det
1325 0 : Int_t info[3] = {0,(Int_t)neg[j].GetNd(),fNlayer[fModule]};
1326 :
1327 0 : lp[2]=0.0085*0.0085; //SigmaY2
1328 0 : lp[3]=1.15*1.15; //SigmaZ2
1329 0 : lp[5]=0.0093;
1330 :
1331 0 : AliITSRecPoint cl2(milab,lp,info);
1332 0 : cl2.SetChargeRatio(1.);
1333 0 : cl2.SetType(80);
1334 0 : clusters.push_back( cl2 );
1335 0 : ncl++;
1336 0 : } // cross is within the detector
1337 : //
1338 : //--------------
1339 :
1340 0 : } // end loop over Pside 1D clusters
1341 :
1342 0 : } // bad Nside module
1343 :
1344 : //---------------------------------------------------------
1345 :
1346 0 : } // use bad channels
1347 :
1348 : //cout<<ncl<<" clusters for this module"<<endl;
1349 :
1350 0 : delete [] cnegative;
1351 0 : delete [] cused1;
1352 0 : delete [] negativepair;
1353 0 : delete [] cpositive;
1354 0 : delete [] cused2;
1355 0 : delete [] positivepair;
1356 :
1357 0 : }
|
