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