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