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