merging RecPoints and ClustersV2. All ClusterFinders produce AliITSRecPoints objects...
[u/mrichter/AliRoot.git] / ITS / AliITSClusterFinderV2SSD.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-2003, 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 //            Implementation of the ITS clusterer V2 class                //
17 //                                                                        //
18 //          Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch            //
19 //                                                                        //
20 ///////////////////////////////////////////////////////////////////////////
21
22
23 #include "AliITSClusterFinderV2SSD.h"
24 #include "AliITSRecPoint.h"
25 #include "AliITSDetTypeRec.h"
26 #include "AliRawReader.h"
27 #include "AliITSRawStreamSSD.h"
28 #include <TClonesArray.h>
29 #include "AliITSdigitSSD.h"
30
31 ClassImp(AliITSClusterFinderV2SSD)
32
33
34 AliITSClusterFinderV2SSD::AliITSClusterFinderV2SSD(AliITSDetTypeRec* dettyp):AliITSClusterFinderV2(dettyp){
35
36   //Default constructor
37
38   fLastSSD1=fDetTypeRec->GetITSgeom()->GetModuleIndex(6,1,1)-1;
39   fYpitchSSD=0.0095;
40   fHwSSD=3.65;
41   fHlSSD=2.00;
42   fTanP=0.0275;
43   fTanN=0.0075;
44
45
46
47 }
48  
49
50 void AliITSClusterFinderV2SSD::FindRawClusters(Int_t mod){
51
52   //Find clusters V2
53   SetModule(mod);
54   FindClustersSSD(fDigits);
55
56 }
57
58 void AliITSClusterFinderV2SSD::FindClustersSSD(TClonesArray *alldigits) {
59   //------------------------------------------------------------
60   // Actual SSD cluster finder
61   //------------------------------------------------------------
62   Int_t smaxall=alldigits->GetEntriesFast();
63   if (smaxall==0) return;
64   TObjArray *digits = new TObjArray;
65   for (Int_t i=0;i<smaxall; i++){
66     AliITSdigitSSD *d=(AliITSdigitSSD*)alldigits->UncheckedAt(i);
67     if (d->GetSignal()<3) continue;
68     digits->AddLast(d);
69   }
70   Int_t smax = digits->GetEntriesFast();
71   if (smax==0) return;
72   
73   const Int_t kMax=1000;
74   Int_t np=0, nn=0; 
75   Ali1Dcluster pos[kMax], neg[kMax];
76   Float_t y=0., q=0., qmax=0.; 
77   Int_t lab[4]={-2,-2,-2,-2};
78   
79   AliITSdigitSSD *d=(AliITSdigitSSD*)digits->UncheckedAt(0);
80   q += d->GetSignal();
81   y += d->GetCoord2()*d->GetSignal();
82   qmax=d->GetSignal();
83   lab[0]=d->GetTrack(0); lab[1]=d->GetTrack(1); lab[2]=d->GetTrack(2);
84   Int_t curr=d->GetCoord2();
85   Int_t flag=d->GetCoord1();
86   Int_t *n=&nn;
87   Ali1Dcluster *c=neg;
88   Int_t nd=1;
89   Int_t milab[10];
90   for (Int_t ilab=0;ilab<10;ilab++){
91     milab[ilab]=-2;
92   }
93   milab[0]=d->GetTrack(0); milab[1]=d->GetTrack(1); milab[2]=d->GetTrack(2);
94
95   for (Int_t s=1; s<smax; s++) {
96       d=(AliITSdigitSSD*)digits->UncheckedAt(s);      
97       Int_t strip=d->GetCoord2();
98       if ((strip-curr) > 1 || flag!=d->GetCoord1()) {
99          c[*n].SetY(y/q);
100          c[*n].SetQ(q);
101          c[*n].SetNd(nd);
102          CheckLabels2(milab);
103          c[*n].SetLabels(milab);
104          //Split suspiciously big cluster
105          /*
106          if (nd>10&&nd<16){
107            c[*n].SetY(y/q-0.3*nd);
108            c[*n].SetQ(0.5*q);
109            (*n)++;
110            if (*n==MAX) {
111              Error("FindClustersSSD","Too many 1D clusters !");
112               return;
113            }
114            c[*n].SetY(y/q-0.0*nd);
115            c[*n].SetQ(0.5*q);
116            c[*n].SetNd(nd);
117            (*n)++;
118            if (*n==MAX) {
119              Error("FindClustersSSD","Too many 1D clusters !");
120               return;
121            }
122            //
123            c[*n].SetY(y/q+0.3*nd);
124            c[*n].SetQ(0.5*q);
125            c[*n].SetNd(nd);
126            c[*n].SetLabels(milab);
127          }
128          else{
129          */
130          if (nd>4&&nd<25) {
131            c[*n].SetY(y/q-0.25*nd);
132            c[*n].SetQ(0.5*q);
133            (*n)++;
134            if (*n==kMax) {
135              Error("FindClustersSSD","Too many 1D clusters !");
136              return;
137            }
138            c[*n].SetY(y/q+0.25*nd);
139            c[*n].SetQ(0.5*q);
140            c[*n].SetNd(nd);
141            c[*n].SetLabels(milab);
142          }       
143          (*n)++;
144          if (*n==kMax) {
145           Error("FindClustersSSD","Too many 1D clusters !");
146           return;
147          }
148          y=q=qmax=0.;
149          nd=0;
150          lab[0]=lab[1]=lab[2]=-2;
151          //
152          for (Int_t ilab=0;ilab<10;ilab++){
153            milab[ilab]=-2;
154          }
155          //
156          if (flag!=d->GetCoord1()) { n=&np; c=pos; }
157       }
158       flag=d->GetCoord1();
159       q += d->GetSignal();
160       y += d->GetCoord2()*d->GetSignal();
161       nd++;
162       if (d->GetSignal()>qmax) {
163          qmax=d->GetSignal();
164          lab[0]=d->GetTrack(0); lab[1]=d->GetTrack(1); lab[2]=d->GetTrack(2);
165       }
166       for (Int_t ilab=0;ilab<10;ilab++) {
167         if (d->GetTrack(ilab)>=0) AddLabel(milab, (d->GetTrack(ilab))); 
168       }
169       curr=strip;
170   }
171   c[*n].SetY(y/q);
172   c[*n].SetQ(q);
173   c[*n].SetNd(nd);
174   c[*n].SetLabels(lab);
175   //Split suspiciously big cluster
176   if (nd>4 && nd<25) {
177      c[*n].SetY(y/q-0.25*nd);
178      c[*n].SetQ(0.5*q);
179      (*n)++;
180      if (*n==kMax) {
181         Error("FindClustersSSD","Too many 1D clusters !");
182         return;
183      }
184      c[*n].SetY(y/q+0.25*nd);
185      c[*n].SetQ(0.5*q);
186      c[*n].SetNd(nd);
187      c[*n].SetLabels(lab);
188   }
189   (*n)++;
190   if (*n==kMax) {
191      Error("FindClustersSSD","Too many 1D clusters !");
192      return;
193   }
194
195   FindClustersSSD(neg, nn, pos, np);
196 }
197
198
199 void AliITSClusterFinderV2SSD::RawdataToClusters(AliRawReader* rawReader,TClonesArray** clusters){
200
201     //------------------------------------------------------------
202   // This function creates ITS clusters from raw data
203   //------------------------------------------------------------
204   rawReader->Reset();
205   AliITSRawStreamSSD inputSSD(rawReader);
206   FindClustersSSD(&inputSSD,clusters);
207   
208 }
209
210 void AliITSClusterFinderV2SSD::FindClustersSSD(AliITSRawStream* input, 
211                                         TClonesArray** clusters) 
212 {
213   //------------------------------------------------------------
214   // Actual SSD cluster finder for raw data
215   //------------------------------------------------------------
216   Int_t nClustersSSD = 0;
217   const Int_t kMax = 1000;
218   Ali1Dcluster clusters1D[2][kMax];
219   Int_t nClusters[2] = {0, 0};
220   Int_t lab[3]={-2,-2,-2};
221   Float_t q = 0.;
222   Float_t y = 0.;
223   Int_t nDigits = 0;
224   Int_t prevStrip = -1;
225   Int_t prevFlag = -1;
226   Int_t prevModule = -1;
227
228   // read raw data input stream
229   while (kTRUE) {
230     Bool_t next = input->Next();
231
232     if(input->GetSignal()<3 && next) continue;
233     // check if a new cluster starts
234     Int_t strip = input->GetCoord2();
235     Int_t flag = input->GetCoord1();
236     if ((!next || (input->GetModuleID() != prevModule)||
237          (strip-prevStrip > 1) || (flag != prevFlag)) &&
238         (nDigits > 0)) {
239       if (nClusters[prevFlag] == kMax) {
240         Error("FindClustersSSD", "Too many 1D clusters !");
241         return;
242       }
243       Ali1Dcluster& cluster = clusters1D[prevFlag][nClusters[prevFlag]++];
244       cluster.SetY(y/q);
245       cluster.SetQ(q);
246       cluster.SetNd(nDigits);
247       cluster.SetLabels(lab);
248
249       //Split suspiciously big cluster
250       if (nDigits > 4&&nDigits < 25) {
251         cluster.SetY(y/q - 0.25*nDigits);
252         cluster.SetQ(0.5*q);
253         if (nClusters[prevFlag] == kMax) {
254           Error("FindClustersSSD", "Too many 1D clusters !");
255           return;
256         }
257         Ali1Dcluster& cluster2 = clusters1D[prevFlag][nClusters[prevFlag]++];
258         cluster2.SetY(y/q + 0.25*nDigits);
259         cluster2.SetQ(0.5*q);
260         cluster2.SetNd(nDigits);
261         cluster2.SetLabels(lab);
262       }
263       y = q = 0.;
264       nDigits = 0;
265     }
266
267     if (!next || (input->GetModuleID() != prevModule)) {
268       Int_t iModule = prevModule;
269
270       // when all data from a module was read, search for clusters
271       if (prevFlag >= 0) {
272         clusters[iModule] = new TClonesArray("AliITSRecPoint");
273         fModule = iModule;
274         FindClustersSSD(&clusters1D[0][0], nClusters[0], 
275                         &clusters1D[1][0], nClusters[1], clusters[iModule]);
276         Int_t nClusters = clusters[iModule]->GetEntriesFast();
277         nClustersSSD += nClusters;
278       }
279
280       if (!next) break;
281       nClusters[0] = nClusters[1] = 0;
282       y = q = 0.;
283       nDigits = 0;
284     }
285
286     // add digit to current cluster
287     q += input->GetSignal();
288     y += strip * input->GetSignal();
289     nDigits++;
290     prevStrip = strip;
291     prevFlag = flag;
292     prevModule = input->GetModuleID();
293
294   }
295
296   Info("FindClustersSSD", "found clusters in ITS SSD: %d", nClustersSSD);
297 }
298
299 void AliITSClusterFinderV2SSD::
300 FindClustersSSD(Ali1Dcluster* neg, Int_t nn, 
301                 Ali1Dcluster* pos, Int_t np,
302                 TClonesArray *clusters) {
303   //------------------------------------------------------------
304   // Actual SSD cluster finder
305   //------------------------------------------------------------
306   TClonesArray &cl=*clusters;
307   //
308   Float_t tanp=fTanP, tann=fTanN;
309   if (fModule>fLastSSD1) {tann=fTanP; tanp=fTanN;}
310   Int_t idet=fNdet[fModule];
311   Int_t ncl=0;
312   //
313   Int_t negativepair[30000];
314   Int_t cnegative[3000];  
315   Int_t cused1[3000];
316   Int_t positivepair[30000];
317   Int_t cpositive[3000];
318   Int_t cused2[3000];
319   for (Int_t i=0;i<3000;i++) {cnegative[i]=0; cused1[i]=0;}
320   for (Int_t i=0;i<3000;i++) {cpositive[i]=0; cused2[i]=0;}
321   static Short_t pairs[1000][1000];
322   memset(pairs,0,sizeof(Short_t)*1000000);
323 //   Short_t ** pairs = new Short_t*[1000];
324 //   for (Int_t i=0; i<1000; i++) {
325 //     pairs[i] = new Short_t[1000];
326 //     memset(pairs[i],0,sizeof(Short_t)*1000);
327 //   }  
328   //
329   // find available pairs
330   //
331   for (Int_t i=0; i<np; i++) {
332     Float_t yp=pos[i].GetY()*fYpitchSSD; 
333     if (pos[i].GetQ()<3) continue;
334     for (Int_t j=0; j<nn; j++) {
335       if (neg[j].GetQ()<3) continue;
336       Float_t yn=neg[j].GetY()*fYpitchSSD;
337       Float_t zt=(2*fHlSSD*tanp + yp - yn)/(tann+tanp);
338       Float_t yt=yn + tann*zt;
339       zt-=fHlSSD; yt-=fHwSSD;
340       if (TMath::Abs(yt)<fHwSSD+0.01)
341       if (TMath::Abs(zt)<fHlSSD+0.01*(neg[j].GetNd()+pos[i].GetNd())) {
342         negativepair[i*10+cnegative[i]] =j;  //index
343         positivepair[j*10+cpositive[j]] =i;
344         cnegative[i]++;  //counters
345         cpositive[j]++; 
346         pairs[i][j]=100;
347       }
348     }
349   }
350   //
351   for (Int_t i=0; i<np; i++) {
352     Float_t yp=pos[i].GetY()*fYpitchSSD; 
353     if (pos[i].GetQ()<3) continue;
354     for (Int_t j=0; j<nn; j++) {
355       if (neg[j].GetQ()<3) continue;
356       if (cpositive[j]&&cnegative[i]) continue;
357       Float_t yn=neg[j].GetY()*fYpitchSSD;
358       Float_t zt=(2*fHlSSD*tanp + yp - yn)/(tann+tanp);
359       Float_t yt=yn + tann*zt;
360       zt-=fHlSSD; yt-=fHwSSD;
361       if (TMath::Abs(yt)<fHwSSD+0.1)
362       if (TMath::Abs(zt)<fHlSSD+0.15) {
363         if (cnegative[i]==0) pos[i].SetNd(100);  // not available pair
364         if (cpositive[j]==0) neg[j].SetNd(100);  // not available pair
365         negativepair[i*10+cnegative[i]] =j;  //index
366         positivepair[j*10+cpositive[j]] =i;
367         cnegative[i]++;  //counters
368         cpositive[j]++; 
369         pairs[i][j]=100;
370       }
371     }
372   }
373   //
374   Float_t lp[5];
375   Int_t milab[10];
376   Double_t ratio;
377   
378   //
379   // sign gold tracks
380   //
381   for (Int_t ip=0;ip<np;ip++){
382     Float_t ybest=1000,zbest=1000,qbest=0;
383     //
384     // select gold clusters
385     if ( (cnegative[ip]==1) && cpositive[negativepair[10*ip]]==1){ 
386       Float_t yp=pos[ip].GetY()*fYpitchSSD; 
387       Int_t j = negativepair[10*ip];      
388       ratio = (pos[ip].GetQ()-neg[j].GetQ())/(pos[ip].GetQ()+neg[j].GetQ());
389       //
390       Float_t yn=neg[j].GetY()*fYpitchSSD;
391       Float_t zt=(2*fHlSSD*tanp + yp - yn)/(tann+tanp);
392       Float_t yt=yn + tann*zt;
393       zt-=fHlSSD; yt-=fHwSSD;
394       ybest=yt; zbest=zt; 
395       qbest=0.5*(pos[ip].GetQ()+neg[j].GetQ());
396       lp[0]=-(-ybest+fYshift[fModule]);
397       lp[1]=  -zbest+fZshift[fModule];
398       lp[2]=0.0025*0.0025;  //SigmaY2
399       lp[3]=0.110*0.110;  //SigmaZ2
400       
401       lp[4]=qbest;        //Q
402       for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2;
403       for (Int_t ilab=0;ilab<3;ilab++){
404         milab[ilab] = pos[ip].GetLabel(ilab);
405         milab[ilab+3] = neg[j].GetLabel(ilab);
406       }
407       //
408       CheckLabels2(milab);
409       milab[3]=(((ip<<10) + j)<<10) + idet; // pos|neg|det
410       Int_t info[3] = {pos[ip].GetNd(),neg[j].GetNd(),fNlayer[fModule]};
411       AliITSRecPoint * cl2;
412       if(clusters){
413         cl2 = new (cl[ncl]) AliITSRecPoint(fModule,fDetTypeRec->GetITSgeom(),milab,lp,info);
414         cl2->SetChargeRatio(ratio);     
415         cl2->SetType(1);
416         pairs[ip][j]=1;
417         if ((pos[ip].GetNd()+neg[j].GetNd())>6){ //multi cluster
418           cl2->SetType(2);
419           pairs[ip][j]=2;
420         }
421         cused1[ip]++;
422         cused2[j]++;
423         
424       }
425       else{
426         cl2 = new AliITSRecPoint(fModule,fDetTypeRec->GetITSgeom(),milab,lp,info);      
427         cl2->SetChargeRatio(ratio);     
428         cl2->SetType(1);
429         pairs[ip][j]=1;
430         if ((pos[ip].GetNd()+neg[j].GetNd())>6){ //multi cluster
431           cl2->SetType(2);
432           pairs[ip][j]=2;
433         }
434         cused1[ip]++;
435         cused2[j]++;
436         fDetTypeRec->AddRecPoint(*cl2);
437       }
438       ncl++;
439     }
440   }
441     
442   for (Int_t ip=0;ip<np;ip++){
443     Float_t ybest=1000,zbest=1000,qbest=0;
444     //
445     //
446     // select "silber" cluster
447     if ( cnegative[ip]==1 && cpositive[negativepair[10*ip]]==2){
448       Int_t in  = negativepair[10*ip];
449       Int_t ip2 = positivepair[10*in];
450       if (ip2==ip) ip2 =  positivepair[10*in+1];
451       Float_t pcharge = pos[ip].GetQ()+pos[ip2].GetQ();
452       if (TMath::Abs(pcharge-neg[in].GetQ())<10){
453         //
454         // add first pair
455         if (pairs[ip][in]==100){  //
456           Float_t yp=pos[ip].GetY()*fYpitchSSD; 
457           Float_t yn=neg[in].GetY()*fYpitchSSD;
458           Float_t zt=(2*fHlSSD*tanp + yp - yn)/(tann+tanp);
459           Float_t yt=yn + tann*zt;
460           zt-=fHlSSD; yt-=fHwSSD;
461           ybest =yt;  zbest=zt; 
462           qbest =pos[ip].GetQ();
463           lp[0]=-(-ybest+fYshift[fModule]);
464           lp[1]=  -zbest+fZshift[fModule];
465           lp[2]=0.0025*0.0025;  //SigmaY2
466           lp[3]=0.110*0.110;  //SigmaZ2
467           
468           lp[4]=qbest;        //Q
469           for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2;
470           for (Int_t ilab=0;ilab<3;ilab++){
471             milab[ilab] = pos[ip].GetLabel(ilab);
472             milab[ilab+3] = neg[in].GetLabel(ilab);
473           }
474           //
475           CheckLabels2(milab);
476           ratio = (pos[ip].GetQ()-neg[in].GetQ())/(pos[ip].GetQ()+neg[in].GetQ());
477           milab[3]=(((ip<<10) + in)<<10) + idet; // pos|neg|det
478           Int_t info[3] = {pos[ip].GetNd(),neg[in].GetNd(),fNlayer[fModule]};
479
480           AliITSRecPoint * cl2;
481           if(clusters){
482             cl2 = new (cl[ncl]) AliITSRecPoint(fModule,fDetTypeRec->GetITSgeom(),milab,lp,info);
483             cl2->SetChargeRatio(ratio);         
484             cl2->SetType(5);
485             pairs[ip][in] = 5;
486             if ((pos[ip].GetNd()+neg[in].GetNd())>6){ //multi cluster
487               cl2->SetType(6);
488               pairs[ip][in] = 6;
489             }       
490           }
491           else{
492             cl2 = new AliITSRecPoint(fModule,fDetTypeRec->GetITSgeom(),milab,lp,info);
493             cl2->SetChargeRatio(ratio);         
494             cl2->SetType(5);
495             pairs[ip][in] = 5;
496             if ((pos[ip].GetNd()+neg[in].GetNd())>6){ //multi cluster
497               cl2->SetType(6);
498               pairs[ip][in] = 6;
499             }
500             
501             fDetTypeRec->AddRecPoint(*cl2);
502           }
503           ncl++;
504         }
505         
506         //
507         // add second pair
508         
509       //        if (!(cused1[ip2] || cused2[in])){  //
510         if (pairs[ip2][in]==100){
511           Float_t yp=pos[ip2].GetY()*fYpitchSSD;
512           Float_t yn=neg[in].GetY()*fYpitchSSD;
513           Float_t zt=(2*fHlSSD*tanp + yp - yn)/(tann+tanp);
514           Float_t yt=yn + tann*zt;
515           zt-=fHlSSD; yt-=fHwSSD;
516           ybest =yt;  zbest=zt; 
517           qbest =pos[ip2].GetQ();
518           lp[0]=-(-ybest+fYshift[fModule]);
519           lp[1]=  -zbest+fZshift[fModule];
520           lp[2]=0.0025*0.0025;  //SigmaY2
521           lp[3]=0.110*0.110;  //SigmaZ2
522           
523           lp[4]=qbest;        //Q
524           for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2;
525           for (Int_t ilab=0;ilab<3;ilab++){
526             milab[ilab] = pos[ip2].GetLabel(ilab);
527             milab[ilab+3] = neg[in].GetLabel(ilab);
528           }
529           //
530           CheckLabels2(milab);
531           ratio = (pos[ip2].GetQ()-neg[in].GetQ())/(pos[ip2].GetQ()+neg[in].GetQ());
532           milab[3]=(((ip2<<10) + in)<<10) + idet; // pos|neg|det
533           Int_t info[3] = {pos[ip2].GetNd(),neg[in].GetNd(),fNlayer[fModule]};
534           
535           AliITSRecPoint * cl2;
536           if(clusters){
537             cl2 = new (cl[ncl]) AliITSRecPoint(fModule,fDetTypeRec->GetITSgeom(),milab,lp,info);
538             cl2->SetChargeRatio(ratio);         
539             cl2->SetType(5);
540             pairs[ip2][in] =5;
541             if ((pos[ip2].GetNd()+neg[in].GetNd())>6){ //multi cluster
542               cl2->SetType(6);
543               pairs[ip2][in] =6;
544             }
545           }
546           else{
547             cl2 = new AliITSRecPoint(fModule,fDetTypeRec->GetITSgeom(),milab,lp,info);
548             cl2->SetChargeRatio(ratio);         
549             cl2->SetType(5);
550             pairs[ip2][in] =5;
551             if ((pos[ip2].GetNd()+neg[in].GetNd())>6){ //multi cluster
552               cl2->SetType(6);
553               pairs[ip2][in] =6;
554             }
555
556             fDetTypeRec->AddRecPoint(*cl2);
557           }
558           ncl++;
559         }       
560         cused1[ip]++;
561         cused1[ip2]++;
562         cused2[in]++;
563       }
564     }    
565   }
566
567   //  
568   for (Int_t jn=0;jn<nn;jn++){
569     if (cused2[jn]) continue;
570     Float_t ybest=1000,zbest=1000,qbest=0;
571     // select "silber" cluster
572     if ( cpositive[jn]==1 && cnegative[positivepair[10*jn]]==2){
573       Int_t ip  = positivepair[10*jn];
574       Int_t jn2 = negativepair[10*ip];
575       if (jn2==jn) jn2 =  negativepair[10*ip+1];
576       Float_t pcharge = neg[jn].GetQ()+neg[jn2].GetQ();
577       //
578       if (TMath::Abs(pcharge-pos[ip].GetQ())<10){
579         //
580         // add first pair
581         //      if (!(cused1[ip]||cused2[jn])){
582         if (pairs[ip][jn]==100){
583           Float_t yn=neg[jn].GetY()*fYpitchSSD; 
584           Float_t yp=pos[ip].GetY()*fYpitchSSD;
585           Float_t zt=(2*fHlSSD*tanp + yp - yn)/(tann+tanp);
586           Float_t yt=yn + tann*zt;
587           zt-=fHlSSD; yt-=fHwSSD;
588           ybest =yt;  zbest=zt; 
589           qbest =neg[jn].GetQ();
590           lp[0]=-(-ybest+fYshift[fModule]);
591           lp[1]=  -zbest+fZshift[fModule];
592           lp[2]=0.0025*0.0025;  //SigmaY2
593           lp[3]=0.110*0.110;  //SigmaZ2
594           
595           lp[4]=qbest;        //Q
596           for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2;
597           for (Int_t ilab=0;ilab<3;ilab++){
598             milab[ilab] = pos[ip].GetLabel(ilab);
599             milab[ilab+3] = neg[jn].GetLabel(ilab);
600           }
601           //
602           CheckLabels2(milab);
603           ratio = (pos[ip].GetQ()-neg[jn].GetQ())/(pos[ip].GetQ()+neg[jn].GetQ());
604           milab[3]=(((ip<<10) + jn)<<10) + idet; // pos|neg|det
605           Int_t info[3] = {pos[ip].GetNd(),neg[jn].GetNd(),fNlayer[fModule]};
606
607           AliITSRecPoint * cl2;
608           if(clusters){
609             cl2 = new (cl[ncl]) AliITSRecPoint(fModule,fDetTypeRec->GetITSgeom(),milab,lp,info);
610             cl2->SetChargeRatio(ratio);         
611             cl2->SetType(7);
612             pairs[ip][jn] =7;
613             if ((pos[ip].GetNd()+neg[jn].GetNd())>6){ //multi cluster
614               cl2->SetType(8);
615               pairs[ip][jn]=8;
616             }
617
618           }
619           else{
620             cl2 = new AliITSRecPoint(fModule,fDetTypeRec->GetITSgeom(),milab,lp,info);
621             cl2->SetChargeRatio(ratio);         
622             cl2->SetType(7);
623             pairs[ip][jn] =7;
624             if ((pos[ip].GetNd()+neg[jn].GetNd())>6){ //multi cluster
625               cl2->SetType(8);
626               pairs[ip][jn]=8;
627             }
628
629             fDetTypeRec->AddRecPoint(*cl2);
630           }
631           ncl++;
632         }
633         //
634         // add second pair
635         //      if (!(cused1[ip]||cused2[jn2])){
636         if (pairs[ip][jn2]==100){
637           Float_t yn=neg[jn2].GetY()*fYpitchSSD; 
638           Double_t yp=pos[ip].GetY()*fYpitchSSD; 
639           Double_t zt=(2*fHlSSD*tanp + yp - yn)/(tann+tanp);
640           Double_t yt=yn + tann*zt;
641           zt-=fHlSSD; yt-=fHwSSD;
642           ybest =yt;  zbest=zt; 
643           qbest =neg[jn2].GetQ();
644           lp[0]=-(-ybest+fYshift[fModule]);
645           lp[1]=  -zbest+fZshift[fModule];
646           lp[2]=0.0025*0.0025;  //SigmaY2
647           lp[3]=0.110*0.110;  //SigmaZ2
648           
649           lp[4]=qbest;        //Q
650           for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2;
651           for (Int_t ilab=0;ilab<3;ilab++){
652             milab[ilab] = pos[ip].GetLabel(ilab);
653             milab[ilab+3] = neg[jn2].GetLabel(ilab);
654           }
655           //
656           CheckLabels2(milab);
657           ratio = (pos[ip].GetQ()-neg[jn2].GetQ())/(pos[ip].GetQ()+neg[jn2].GetQ());
658           milab[3]=(((ip<<10) + jn2)<<10) + idet; // pos|neg|det
659           Int_t info[3] = {pos[ip].GetNd(),neg[jn2].GetNd(),fNlayer[fModule]};
660           AliITSRecPoint * cl2;
661           if(clusters){
662             cl2 = new (cl[ncl]) AliITSRecPoint(fModule,fDetTypeRec->GetITSgeom(),milab,lp,info);
663             cl2->SetChargeRatio(ratio);         
664             pairs[ip][jn2]=7;
665             cl2->SetType(7);
666             if ((pos[ip].GetNd()+neg[jn2].GetNd())>6){ //multi cluster
667               cl2->SetType(8);
668               pairs[ip][jn2]=8;
669             }
670             
671           }
672           else{
673             cl2 = new AliITSRecPoint(fModule,fDetTypeRec->GetITSgeom(),milab,lp,info);
674             cl2->SetChargeRatio(ratio);         
675             pairs[ip][jn2]=7;
676             cl2->SetType(7);
677             if ((pos[ip].GetNd()+neg[jn2].GetNd())>6){ //multi cluster
678               cl2->SetType(8);
679               pairs[ip][jn2]=8;
680             }
681             
682             fDetTypeRec->AddRecPoint(*cl2);
683           }
684
685           ncl++;
686         }
687         cused1[ip]++;
688         cused2[jn]++;
689         cused2[jn2]++;
690       }
691     }    
692   }
693   
694   for (Int_t ip=0;ip<np;ip++){
695     Float_t ybest=1000,zbest=1000,qbest=0;
696     //
697     // 2x2 clusters
698     //
699     if ( (cnegative[ip]<5) && cpositive[negativepair[10*ip]]<5){ 
700       Float_t minchargediff =4.;
701       Int_t j=-1;
702       for (Int_t di=0;di<cnegative[ip];di++){
703         Int_t   jc = negativepair[ip*10+di];
704         Float_t chargedif = pos[ip].GetQ()-neg[jc].GetQ();
705         if (TMath::Abs(chargedif)<minchargediff){
706           j =jc;
707           minchargediff = TMath::Abs(chargedif);
708         }
709       }
710       if (j<0) continue;  // not proper cluster      
711       Int_t count =0;
712       for (Int_t di=0;di<cnegative[ip];di++){
713         Int_t   jc = negativepair[ip*10+di];
714         Float_t chargedif = pos[ip].GetQ()-neg[jc].GetQ();
715         if (TMath::Abs(chargedif)<minchargediff+3.) count++;
716       }
717       if (count>1) continue;  // more than one "proper" cluster for positive
718       //
719       count =0;
720       for (Int_t dj=0;dj<cpositive[j];dj++){
721         Int_t   ic  = positivepair[j*10+dj];
722         Float_t chargedif = pos[ic].GetQ()-neg[j].GetQ();
723         if (TMath::Abs(chargedif)<minchargediff+3.) count++;
724       }
725       if (count>1) continue;  // more than one "proper" cluster for negative
726       
727       Int_t jp = 0;
728       
729       count =0;
730       for (Int_t dj=0;dj<cnegative[jp];dj++){
731         Int_t   ic = positivepair[jp*10+dj];
732         Float_t chargedif = pos[ic].GetQ()-neg[jp].GetQ();
733         if (TMath::Abs(chargedif)<minchargediff+4.) count++;
734       }
735       if (count>1) continue;   
736       if (pairs[ip][j]<100) continue;
737       //
738       //almost gold clusters
739       Float_t yp=pos[ip].GetY()*fYpitchSSD; 
740       Float_t yn=neg[j].GetY()*fYpitchSSD;
741       Float_t zt=(2*fHlSSD*tanp + yp - yn)/(tann+tanp);
742       Float_t yt=yn + tann*zt;
743       zt-=fHlSSD; yt-=fHwSSD;
744       ybest=yt; zbest=zt; 
745       qbest=0.5*(pos[ip].GetQ()+neg[j].GetQ());
746       lp[0]=-(-ybest+fYshift[fModule]);
747       lp[1]=  -zbest+fZshift[fModule];
748       lp[2]=0.0025*0.0025;  //SigmaY2
749       lp[3]=0.110*0.110;  //SigmaZ2     
750       lp[4]=qbest;        //Q
751       for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2;
752       for (Int_t ilab=0;ilab<3;ilab++){
753         milab[ilab] = pos[ip].GetLabel(ilab);
754         milab[ilab+3] = neg[j].GetLabel(ilab);
755       }
756       //
757       CheckLabels2(milab);
758       ratio = (pos[ip].GetQ()-neg[j].GetQ())/(pos[ip].GetQ()+neg[j].GetQ());
759       milab[3]=(((ip<<10) + j)<<10) + idet; // pos|neg|det
760       Int_t info[3] = {pos[ip].GetNd(),neg[j].GetNd(),fNlayer[fModule]};
761       AliITSRecPoint * cl2;
762       if(clusters){
763         cl2 = new (cl[ncl]) AliITSRecPoint(fModule,fDetTypeRec->GetITSgeom(),milab,lp,info);
764         cl2->SetChargeRatio(ratio);     
765         cl2->SetType(10);
766         pairs[ip][j]=10;
767         if ((pos[ip].GetNd()+neg[j].GetNd())>6){ //multi cluster
768           cl2->SetType(11);
769           pairs[ip][j]=11;
770         }
771         cused1[ip]++;
772         cused2[j]++;      
773       }
774       else{
775         cl2 = new AliITSRecPoint(fModule,fDetTypeRec->GetITSgeom(),milab,lp,info);
776         cl2->SetChargeRatio(ratio);     
777         cl2->SetType(10);
778         pairs[ip][j]=10;
779         if ((pos[ip].GetNd()+neg[j].GetNd())>6){ //multi cluster
780           cl2->SetType(11);
781           pairs[ip][j]=11;
782         }
783         cused1[ip]++;
784         cused2[j]++;      
785         
786         fDetTypeRec->AddRecPoint(*cl2);
787       }      
788       ncl++;
789     }
790
791   }
792   
793   //  
794   for (Int_t i=0; i<np; i++) {
795     Float_t ybest=1000,zbest=1000,qbest=0;
796     Float_t yp=pos[i].GetY()*fYpitchSSD; 
797     if (pos[i].GetQ()<3) continue;
798     for (Int_t j=0; j<nn; j++) {
799     //    for (Int_t di = 0;di<cpositive[i];di++){
800     //  Int_t j = negativepair[10*i+di];
801       if (neg[j].GetQ()<3) continue;
802       if (cused2[j]||cused1[i]) continue;      
803       if (pairs[i][j]>0 &&pairs[i][j]<100) continue;
804       ratio = (pos[i].GetQ()-neg[j].GetQ())/(pos[i].GetQ()+neg[j].GetQ());      
805       Float_t yn=neg[j].GetY()*fYpitchSSD;
806       Float_t zt=(2*fHlSSD*tanp + yp - yn)/(tann+tanp);
807       Float_t yt=yn + tann*zt;
808       zt-=fHlSSD; yt-=fHwSSD;
809       if (TMath::Abs(yt)<fHwSSD+0.01)
810       if (TMath::Abs(zt)<fHlSSD+0.01*(neg[j].GetNd()+pos[i].GetNd())) {
811         ybest=yt; zbest=zt; 
812         qbest=0.5*(pos[i].GetQ()+neg[j].GetQ());
813         lp[0]=-(-ybest+fYshift[fModule]);
814         lp[1]=  -zbest+fZshift[fModule];
815         lp[2]=0.0025*0.0025;  //SigmaY2
816         lp[3]=0.110*0.110;  //SigmaZ2
817
818         lp[4]=qbest;        //Q
819         for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2;
820         for (Int_t ilab=0;ilab<3;ilab++){
821           milab[ilab] = pos[i].GetLabel(ilab);
822           milab[ilab+3] = neg[j].GetLabel(ilab);
823         }
824         //
825         CheckLabels2(milab);
826         milab[3]=(((i<<10) + j)<<10) + idet; // pos|neg|det
827         Int_t info[3] = {pos[i].GetNd(),neg[j].GetNd(),fNlayer[fModule]};
828         AliITSRecPoint * cl2;
829         if(clusters){
830           cl2 = new (cl[ncl]) AliITSRecPoint(fModule,fDetTypeRec->GetITSgeom(),milab,lp,info);
831           cl2->SetChargeRatio(ratio);
832           cl2->SetType(100+cpositive[j]+cnegative[i]);    
833         }
834         else{
835           cl2 = new AliITSRecPoint(fModule,fDetTypeRec->GetITSgeom(),milab,lp,info);
836           cl2->SetChargeRatio(ratio);
837           cl2->SetType(100+cpositive[j]+cnegative[i]);
838           fDetTypeRec->AddRecPoint(*cl2);
839         }
840         ncl++;
841         //cl2->SetType(0);
842         /*
843           if (pairs[i][j]<100){
844           printf("problem:- %d\n", pairs[i][j]);
845           }
846           if (cnegative[i]<2&&cpositive[j]<2){
847           printf("problem:- %d\n", pairs[i][j]);
848           }
849         */
850       }
851     }
852   }
853
854 //   for (Int_t i=0; i<1000; i++) delete [] pairs[i];
855 //   delete [] pairs;
856
857 }
858
859
860