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