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