c3b894e6b5deea5d103e8ee652cb50979e74fac7
[u/mrichter/AliRoot.git] / ITS / AliITSClusterFinderSSD.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-1999, 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 /*
17
18 Adding rekonstruction facilities
19 Piotr Krzysztof Skowronski 
20 December 1999.
21 */
22
23 /*
24 15 -18 V 2000
25 Eroor counting routines
26 Automatic combination routines improved (traps)
27
28 */
29
30 #include "AliRun.h"
31 #include "AliITSClusterFinderSSD.h"
32  
33
34 const Int_t debug=0;
35
36 ClassImp(AliITSClusterFinderSSD)
37
38 //____________________________________________________________________
39 //
40 //  Constructor
41 //____________________________________________________________________
42 //                                     
43
44
45 AliITSClusterFinderSSD::AliITSClusterFinderSSD(AliITSsegmentation *seg, TClonesArray *digits, TClonesArray *recp)   
46 {
47
48     fSegmentation=seg;
49     fDigits=digits;
50     fRecPoints=recp;
51     
52     fITS=(AliITS*)gAlice->GetModule("ITS");
53     
54     fClusterP =  new TClonesArray ("AliITSclusterSSD",200);    
55     fNClusterP =0;   
56     
57     fClusterN=  new TClonesArray ("AliITSclusterSSD",200);   
58     fNClusterN =0; 
59
60     fPackages =  new TClonesArray ("AliITSpackageSSD",200);    //packages  
61     fNPackages = 0;
62
63         
64     fDigitsIndexP      =  new TArrayI(300);
65     fNDigitsP          =  0;
66     
67     fDigitsIndexN      =  new TArrayI(300);
68     fNDigitsN          =  0;
69     
70     SetAlpha1(1000);
71     SetAlpha2(1000);
72     SetAlpha3(1000);
73
74     
75     fPitch = fSegmentation->Dpx(0);
76     Float_t StereoP,StereoN;
77     fSegmentation->Angles(StereoP,StereoN);
78     fTanP=TMath::Tan(StereoP);
79     fTanN=TMath::Tan(StereoN);
80     fPNsignalRatio=7./8.;         // warning: hard-wired number
81 }
82
83 //-------------------------------------------------------
84 AliITSClusterFinderSSD::~AliITSClusterFinderSSD() {
85    
86     delete fClusterP;
87     delete fClusterN;        
88     delete fPackages;        
89     delete fDigitsIndexP;        
90     delete fDigitsIndexN; 
91 }
92
93 //-------------------------------------------------------
94 void AliITSClusterFinderSSD::InitReconstruction()
95 {
96
97   register Int_t i; //iterator
98
99   for(i=0;i<fNClusterP;i++)
100     {
101       fClusterP->RemoveAt(i);
102     }
103   fNClusterP  =0;
104   for(i=0;i<fNClusterN;i++)
105     {
106       fClusterN->RemoveAt(i);
107     }
108   fNClusterN=0;
109
110   for(i=0;i<fNPackages;i++)
111     {
112       fPackages->RemoveAt(i);
113     }
114
115   fNPackages = 0;
116   fNDigitsP=0;
117   fNDigitsN=0;
118
119   Float_t StereoP,StereoN;
120   fSegmentation->Angles(StereoP,StereoN);
121
122   CalcStepFactor (StereoP,StereoN);
123
124   if (debug) cout<<"fSFF = "<<fSFF<<"  fSFB = "<<fSFB<<"\n";
125 }
126
127
128 //---------------------------------------------
129 void AliITSClusterFinderSSD::FindRawClusters() 
130 {
131
132
133
134 //Piotr Krzysztof Skowronski
135 //Warsaw University of Technology
136 //skowron@if.pw.edu.pl
137
138 // This function findes out all clusters belonging to one module
139 // 1. Zeroes all space after previous module reconstruction
140 // 2. Finds all neighbouring digits
141 // 3. If necesery, resolves for each group of neighbouring digits 
142 //    how many clusters creates it.
143 // 4. Creates packages  
144 // 5. Creates clusters 
145
146   InitReconstruction();  //ad. 1
147   FillDigitsIndex();
148   SortDigits();
149   FindNeighbouringDigits(); //ad. 2
150   SeparateOverlappedClusters();  //ad. 3
151   ClustersToPackages();  //ad. 4
152   ConsumeClusters();
153   PackagesToPoints();   //ad. 5
154   ReconstructNotConsumedClusters();
155  
156 }
157
158
159 //-------------------------------------------------
160 void AliITSClusterFinderSSD::FindNeighbouringDigits()
161 {
162
163
164 //Piotr Krzysztof Skowronski
165 //Warsaw University of Technology
166 //skowron@if.pw.edu.pl
167
168  register Int_t i;
169
170     //If there are any digits on this side, create 1st Cluster,
171     // add to it this digit, and increment number of clusters
172
173  if ((fNDigitsP>0 )  && (fNDigitsN > 0 )) {     
174
175    Int_t currentstripNo;
176    Int_t *dbuffer = new Int_t [300];   //buffer for strip numbers
177    Int_t dnumber;    //curent number of digits in buffer
178    TArrayI &lDigitsIndexP = *fDigitsIndexP;
179    TArrayI &lDigitsIndexN = *fDigitsIndexN;
180    TObjArray &lDigits=*fDigits;
181    TClonesArray &lClusterP = *fClusterP;
182    TClonesArray &lClusterN = *fClusterN;
183   
184    //process P side 
185    dnumber = 1;
186    dbuffer[0]=lDigitsIndexP[0];
187    //If next digit is a neighbour of previous, adds to last cluster this digit
188    for(i=1; i<fNDigitsP; i++) {
189      //reads new digit
190      currentstripNo = ((AliITSdigitSSD*)lDigits[lDigitsIndexP[i]])->
191                                                             GetStripNumber(); 
192      //check if it is a neighbour of a previous one
193      if ( (((AliITSdigitSSD*)lDigits[lDigitsIndexP[i-1]])->GetStripNumber()) 
194            ==  (currentstripNo-1) ) dbuffer[dnumber++]=lDigitsIndexP[i];
195      else  { 
196        //create a new one side cluster
197        new(lClusterP[fNClusterP++]) AliITSclusterSSD(dnumber,dbuffer,fDigits,SIDEP); 
198        dbuffer[0]=lDigitsIndexP[i];
199        dnumber = 1;
200      }
201    } // end loop over fNDigitsP
202    new(lClusterP[fNClusterP++]) AliITSclusterSSD(dnumber,dbuffer,fDigits,SIDEP);
203
204
205    //process N side 
206    //for comments, see above
207    dnumber = 1;
208    dbuffer[0]=lDigitsIndexN[0];
209    //If next digit is a neighbour of previous, adds to last cluster this digit
210    for(i=1; i<fNDigitsN; i++) { 
211      currentstripNo = ((AliITSdigitSSD*)(lDigits[lDigitsIndexN[i]]))->
212                                                             GetStripNumber();
213      if ( (((AliITSdigitSSD*)lDigits[lDigitsIndexN[i-1]])->GetStripNumber()) 
214            == (currentstripNo-1) ) dbuffer[dnumber++]=lDigitsIndexN[i];
215      else {
216        new(lClusterN[fNClusterN++]) AliITSclusterSSD(dnumber,dbuffer,fDigits,SIDEN);
217        dbuffer[0]=lDigitsIndexN[i];
218        dnumber = 1;
219      }
220    } // end loop over fNDigitsN
221    new(lClusterN[fNClusterN++]) AliITSclusterSSD(dnumber,dbuffer,fDigits,SIDEN);
222    delete [] dbuffer; 
223  
224  } // end condition on  NDigits 
225
226  if (debug) cout<<"\n Found clusters: fNClusterP = "<<fNClusterP<<"  fNClusterN ="<<fNClusterN<<"\n";
227
228 }  
229 /**********************************************************************/
230
231
232 void AliITSClusterFinderSSD::SeparateOverlappedClusters()
233 {
234 //************************************************
235 //Piotr Krzysztof Skowronski
236 //Warsaw University of Technology
237 //skowron@if.pw.edu.pl
238
239   register Int_t i,j; //iterator
240
241   Float_t  factor=0.75;            // How many percent must be lower signal 
242                                    // on the middle one digit
243                                    // from its neighbours
244   Int_t    signal0;                //signal on the strip before the current one
245   Int_t    signal1;                //signal on the current one signal
246   Int_t    signal2;                //signal on the strip after the current one
247   TArrayI *splitlist;              //  List of splits
248   Int_t    numerofsplits=0;        // number of splits
249   Int_t    initPsize = fNClusterP; //initial size of the arrays 
250   Int_t    initNsize = fNClusterN; //we have to keep it because it will grow 
251                                    // in this function and it doasn't make 
252                                    // sense to pass through it again
253
254   splitlist = new TArrayI(300);
255
256   for(i=0;i<initPsize;i++)
257   {
258     if (( ((AliITSclusterSSD*)(*fClusterP)[i])->GetNumOfDigits())==1) continue;
259     if (( ((AliITSclusterSSD*)(*fClusterP)[i])->GetNumOfDigits())==2) continue;
260         Int_t nj=(((AliITSclusterSSD*)(*fClusterP)[i])->GetNumOfDigits()-1);
261         for(j=1; j<nj; j++)
262           {
263             signal1=((AliITSclusterSSD*)(*fClusterP)[i])->GetDigitSignal(j);
264             signal0=((AliITSclusterSSD*)(*fClusterP)[i])->GetDigitSignal(j-1);
265             signal2=((AliITSclusterSSD*)(*fClusterP)[i])->GetDigitSignal(j+1);
266             //if signal is less then factor*signal of its neighbours
267             if (  (signal1<(factor*signal0)) && (signal1<(factor*signal2)) )
268              {                                                               
269                (*splitlist)[numerofsplits++]=j;  
270              }
271           } // end loop over number of digits
272           //split this cluster if necessary
273           if(numerofsplits>0) SplitCluster(splitlist,numerofsplits,i,SIDEP); 
274           numerofsplits=0;
275
276           //in signed places (splitlist)
277   } // end loop over clusters on Pside
278
279   for(i=0;i<initNsize;i++) {
280     if (( ((AliITSclusterSSD*)(*fClusterN)[i])->GetNumOfDigits())==1) continue;
281     if (( ((AliITSclusterSSD*)(*fClusterN)[i])->GetNumOfDigits())==2) continue;
282        Int_t nj=(((AliITSclusterSSD*)(*fClusterN)[i])->GetNumOfDigits()-1);
283        for(j=1; j<nj; j++)
284           {
285             signal1=((AliITSclusterSSD*)(*fClusterN)[i])->GetDigitSignal(j);
286             signal0=((AliITSclusterSSD*)(*fClusterN)[i])->GetDigitSignal(j-1);
287             signal2=((AliITSclusterSSD*)(*fClusterN)[i])->GetDigitSignal(j+1);
288             //if signal is less then factor*signal of its neighbours
289             if (  (signal1<(factor*signal0)) && (signal1<(factor*signal2)) ) 
290                (*splitlist)[numerofsplits++]=j;  
291           } // end loop over number of digits 
292           //split this cluster into more clusters
293           if(numerofsplits>0) SplitCluster(splitlist,numerofsplits,i,SIDEN);
294           numerofsplits=0;                                                              //in signed places (splitlist)
295   } // end loop over clusters on Nside
296
297   delete splitlist;
298 }
299
300 //-------------------------------------------------------
301 void AliITSClusterFinderSSD::SplitCluster(TArrayI *list, Int_t nsplits, Int_t index, Bool_t side)
302 {
303
304
305 //Piotr Krzysztof Skowronski
306 //Warsaw University of Technology
307 //skowron@if.pw.edu.pl
308
309   //This function splits one side cluster into more clusters
310   //number of splits is defined by "nsplits"
311   //Place of splits are defined in the TArray "list" 
312   
313   // For further optimisation: Replace this function by two 
314   // specialised ones (each for one side)
315   // save one "if"
316
317   //For comlete comments see AliITSclusterSSD::SplitCluster
318
319
320   register Int_t i; //iterator
321
322   AliITSclusterSSD* curentcluster;
323   Int_t   *tmpdigits = new Int_t[100];
324   Int_t    NN;
325   // side true means P side
326   if (side) {
327      curentcluster =((AliITSclusterSSD*)((*fClusterP)[index])) ;
328      for(i = nsplits; i>0 ;i--) {  
329          NN=curentcluster->SplitCluster((*list)[(i-1)],tmpdigits);
330          new ((*fClusterP)[fNClusterP]) AliITSclusterSSD(NN,tmpdigits,fDigits,side);
331          ( (AliITSclusterSSD*)((*fClusterP)[fNClusterP]) )->
332                                                       SetLeftNeighbour(kTRUE);
333          //if left cluster had neighbour on the right before split 
334          //new should have it too
335          if ( curentcluster->GetRightNeighbour() ) 
336                      ( (AliITSclusterSSD*)((*fClusterP)[fNClusterP]) )->
337                                                      SetRightNeighbour(kTRUE);
338          else curentcluster->SetRightNeighbour(kTRUE); 
339          fNClusterP++;
340      } // end loop over nplits
341   } else {
342      curentcluster =((AliITSclusterSSD*)((*fClusterN)[index]));
343      for(i = nsplits; i>0 ;i--) {  
344          NN=curentcluster->SplitCluster((*list)[(i-1)],tmpdigits);
345          new ((*fClusterN)[fNClusterN]) AliITSclusterSSD(NN,tmpdigits,fDigits,side);
346          ((AliITSclusterSSD*)((*fClusterN)[fNClusterN]))->
347                                                     SetRightNeighbour(kTRUE);
348          if (curentcluster->GetRightNeighbour())
349                       ( (AliITSclusterSSD*)( (*fClusterN)[fNClusterN]) )->
350                                                      SetRightNeighbour(kTRUE);
351          else curentcluster->SetRightNeighbour(kTRUE);      
352          fNClusterN++;
353      } // end loop over nplits
354   } // end if side
355   delete []tmpdigits;
356
357 }
358
359
360 //-------------------------------------------------
361 Int_t AliITSClusterFinderSSD::SortDigitsP(Int_t start, Int_t end)
362 {
363
364
365 //Piotr Krzysztof Skowronski
366 //Warsaw University of Technology
367 //skowron@if.pw.edu.pl
368   
369   Int_t right;
370   Int_t left;
371   if (start != (end - 1) ) 
372     {
373       left=this->SortDigitsP(start,(start+end)/2);
374       right=this->SortDigitsP((start+end)/2,end);  
375       return (left || right);
376     }
377   else
378    { 
379     left =  ((AliITSdigitSSD*)((*fDigits)[(*fDigitsIndexP)[start]]))->GetStripNumber();
380     right= ((AliITSdigitSSD*)((*fDigits)[(*fDigitsIndexP)[end]]))->GetStripNumber();  
381     if( left > right )
382      {
383        Int_t tmp = (*fDigitsIndexP)[start];
384        (*fDigitsIndexP)[start]=(*fDigitsIndexP)[end];
385        (*fDigitsIndexP)[end]=tmp;
386        return 1;
387      }
388     else return 0;
389    }
390 }
391
392
393 //--------------------------------------------------
394
395 Int_t AliITSClusterFinderSSD::SortDigitsN(Int_t start, Int_t end)
396 {
397
398
399 //Piotr Krzysztof Skowronski
400 //Warsaw University of Technology
401 //skowron@if.pw.edu.pl
402
403   Int_t right;
404   Int_t left;
405   if (start != (end - 1)) 
406     {
407       left=this->SortDigitsN(start,(start+end)/2);
408       right=this->SortDigitsN((start+end)/2,end);  
409       return (left || right);
410     }
411   else 
412    {
413     left =((AliITSdigitSSD*)((*fDigits)[(*fDigitsIndexN)[start]]))->GetStripNumber(); 
414     right=((AliITSdigitSSD*)((*fDigits)[(*fDigitsIndexN)[end]]))->GetStripNumber();
415     if ( left > right )
416       {
417         Int_t tmp = (*fDigitsIndexN)[start];
418         (*fDigitsIndexN)[start]=(*fDigitsIndexN)[end];
419         (*fDigitsIndexN)[end]=tmp;
420         return 1;
421       }else return 0;
422    }  
423 }
424
425
426 //------------------------------------------------
427 void AliITSClusterFinderSSD::FillDigitsIndex()
428 {
429
430
431 //Piotr Krzysztof Skowronski
432 //Warsaw University of Technology
433 //skowron@if.pw.edu.pl
434
435  //Fill the indexes of the clusters belonging to a given ITS module
436  //Created by Piotr K. Skowronski, August 7 1999
437
438  Int_t PNs=0, NNs=0;
439  Int_t tmp,bit,k;
440  Int_t N;
441  Int_t i;
442  
443  N = fDigits->GetEntriesFast();
444
445  Int_t* PSidx = new Int_t [N*sizeof(Int_t)];
446  Int_t* NSidx = new Int_t [N*sizeof(Int_t)]; 
447  if (fDigitsIndexP==NULL) fDigitsIndexP = new TArrayI(N);
448  if (fDigitsIndexN==NULL) fDigitsIndexN = new TArrayI(N);
449  
450  AliITSdigitSSD *dig; 
451  
452  for(i = 0 ; i< N; i++ ) {
453       dig=(AliITSdigitSSD*)fDigits->UncheckedAt(i);
454       if(dig->IsSideP()) { 
455            bit=1;
456            tmp=dig->GetStripNumber();
457            // I find this totally unnecessary - it's just a 
458            // CPU consuming double check
459            for(k=0;k<PNs;k++)
460             {
461              if (tmp==PSidx[k])
462               {
463                  if (debug) cout<<"Such a digit exists \n";
464                  bit=0;
465               }
466            }   
467            // end comment 
468         if(bit) {
469             fDigitsIndexP->AddAt(i,fNDigitsP++);
470             PSidx[PNs++]=tmp;
471         }
472       } else {
473          bit=1;
474          tmp=dig->GetStripNumber();
475          // same as above
476          for(k=0;k<NNs;k++)
477           {
478            if (tmp==NSidx[k])
479             {
480              if (debug) cout<<"Such a digit exists \n";
481              bit=0;
482             }
483           } 
484          // end comment
485          if (bit) {
486           fDigitsIndexN->AddAt(i,fNDigitsN++);
487           NSidx[NNs++] =tmp;
488          }
489       }
490  }
491    
492
493  if (debug) cout<<"Digits :  P = "<<fNDigitsP<<"   N = "<<fNDigitsN<<endl;
494
495 }
496
497
498 //-------------------------------------------
499
500 void AliITSClusterFinderSSD::SortDigits()
501 {
502
503
504 //Piotr Krzysztof Skowronski
505 //Warsaw University of Technology
506 //skowron@if.pw.edu.pl
507
508
509
510   Int_t i;
511   if(fNDigitsP>1) 
512   for(i=0;i<fNDigitsP-1;i++)
513   if (SortDigitsP(0,(fNDigitsP-1-i))==0) break;
514
515   if(fNDigitsN>1) 
516     for(i=0;i<fNDigitsN-1;i++)
517   if(SortDigitsN(0,(fNDigitsN-1-i))==0) break;
518 }
519
520
521
522 //----------------------------------------------
523 void AliITSClusterFinderSSD::FillClIndexArrays(Int_t* arrayP, Int_t *arrayN)
524 {
525
526
527 //Piotr Krzysztof Skowronski
528 //Warsaw University of Technology
529 //skowron@if.pw.edu.pl
530
531
532   register Int_t i;
533   for(i=0; i<fNClusterP;i++)
534     {
535       arrayP[i]=i;
536     }
537   for(i=0; i<fNClusterN;i++)
538     {
539       arrayN[i]=i;
540     }
541 }
542
543
544 //------------------------------------------------------
545 void AliITSClusterFinderSSD::SortClusters(Int_t* arrayP, Int_t *arrayN)
546 {
547
548
549 //Piotr Krzysztof Skowronski
550 //Warsaw University of Technology
551 //skowron@if.pw.edu.pl
552
553
554   Int_t i;
555   if(fNClusterP>1) 
556     for(i=0;i<fNClusterP-1;i++)
557       if (SortClustersP(0,(fNClusterP-1),arrayP)==0)  break;
558     
559     
560   if(fNClusterN>1) 
561     for(i=0;i<fNClusterN-1;i++)
562       if (SortClustersN(0,(fNClusterN-1),arrayN)==0)  break;
563
564 }
565
566
567
568 //---------------------------------------------------
569 Int_t AliITSClusterFinderSSD::SortClustersP(Int_t start, Int_t end, Int_t *array)
570 {
571
572
573 //Piotr Krzysztof Skowronski
574 //Warsaw University of Technology
575 //skowron@if.pw.edu.pl
576
577
578
579   Int_t right;
580   Int_t left;
581   if (start != (end - 1) ) {
582       left=this->SortClustersP(start,(start+end)/2,array);
583       right=this->SortClustersP((start+end)/2,end,array);  
584       return (left || right);
585   } else {
586       left =((AliITSclusterSSD*)((*fClusterP)[array[start]]))->
587                                                          GetDigitStripNo(0);
588       right=((AliITSclusterSSD*)((*fClusterP)[array[ end ]]))->
589                                                          GetDigitStripNo(0);
590       if(left>right) {
591          Int_t tmp = array[start];
592          array[start]=array[end];
593          array[end]=tmp;
594          return 1;
595       } else return 0;
596   }
597
598  
599 }
600
601
602
603 //-------------------------------------------------------
604 Int_t AliITSClusterFinderSSD::SortClustersN(Int_t start, Int_t end, Int_t *array)
605 {
606
607
608 //Piotr Krzysztof Skowronski
609 //Warsaw University of Technology
610 //skowron@if.pw.edu.pl
611
612   Int_t right;
613   Int_t left;
614   
615   if (start != (end - 1) ) {
616       left=this->SortClustersN(start,(start+end)/2,array);
617       right=this->SortClustersN((start+end)/2,end,array);  
618       return (left || right);
619   } else {
620       left =((AliITSclusterSSD*)((*fClusterN)[array[start]]))->
621                                                          GetDigitStripNo(0);
622       right=((AliITSclusterSSD*)((*fClusterN)[array[ end ]]))->
623                                                          GetDigitStripNo(0);
624       if( left > right) {
625          Int_t tmp = array[start];
626          array[start]=array[end];
627          array[end]=tmp;
628          return 1;
629       } else return 0;
630     } 
631
632 }
633     
634
635
636 //-------------------------------------------------------
637 void AliITSClusterFinderSSD::ClustersToPackages()
638 {  
639
640
641
642 //Piotr Krzysztof Skowronski
643 //Warsaw University of Technology
644 //skowron@if.pw.edu.pl
645   
646     
647   Int_t *oneSclP = new Int_t[fNClusterP]; //I want to have sorted 1 S clusters
648   Int_t *oneSclN = new Int_t[fNClusterN]; //I can not sort it in TClonesArray
649                                           //so, I create table of indexes and 
650                                           //sort it
651                                           //I do not use TArrayI on purpose
652                                           // MB: well, that's not true that one
653                                           //cannot sort objects in TClonesArray
654   
655   Int_t i,k;    //iterator
656   AliITSpackageSSD *currentpkg;
657   AliITSclusterSSD *currentP;
658   AliITSclusterSSD *currentN;
659   AliITSclusterSSD *fcurrN;
660
661   Int_t     lastNclIndex=0;          //index of last N sidecluster
662   Int_t     Nclpack=0;               //number of P clusters in current package
663   Int_t     Pclpack=0;               //number of N clusters in current package
664   Int_t     firstStripP;
665   Int_t     lastStripP;
666   Int_t     firstStripN;
667   Int_t     lastStripN;
668   Int_t     tmplastNclIndex;
669
670 //Fills in One Side Clusters Index Array
671   FillClIndexArrays(oneSclP,oneSclN); 
672 //Sorts filled Arrays
673   SortClusters(oneSclP,oneSclN);                   
674
675
676   fNPackages=1;      
677   new ((*fPackages)[0]) AliITSpackageSSD(fClusterP,fClusterN);
678   currentpkg = (AliITSpackageSSD*)((*fPackages)[0]);
679   
680   //main loop on sorted P side clusters 
681   for(i=0;i<fNClusterP;i++) {  
682       //Take new P side cluster
683       currentP = GetPSideCluster(oneSclP[i]);
684       currentN = GetNSideCluster(oneSclN[lastNclIndex]);
685       // take a new P cluster or not ?
686       if(IsCrossing(currentP,currentN)) {
687           // don't take a new P cluster
688           Pclpack++;
689           currentP->AddCross(oneSclN[lastNclIndex]);
690           currentN->AddCross(oneSclP[i]);
691           currentpkg->AddPSideCluster(oneSclP[i]);
692           if (Nclpack==0) {
693               currentpkg->AddNSideCluster(oneSclN[lastNclIndex]);          
694               Nclpack++;
695           }  
696           //check if previous N side clusters crosses with it too    
697           for(k=1;k<fSFB+1;k++) {
698               //check if we are not out of array 
699               if ((lastNclIndex-k)>-1) {
700                   fcurrN = GetNSideCluster( oneSclN[lastNclIndex-k] );
701                   if( IsCrossing(currentP,fcurrN) ) {
702                       currentP->AddCross(oneSclN[lastNclIndex-k]);
703                       fcurrN->AddCross(oneSclP[i]);      
704                  } else break; //There is no sense to check rest of clusters
705               } else break;
706           }
707           tmplastNclIndex=lastNclIndex;  
708           //Check if next N side clusters crosses with it too
709           for(k=1;k<fSFF+1;k++) {
710               //check if we are not out of array 
711               if ((tmplastNclIndex+k)<fNClusterN) {
712                   fcurrN = GetNSideCluster( oneSclN[tmplastNclIndex+k] );
713                   if(IsCrossing(currentP,fcurrN) ) {
714                      lastNclIndex++;
715                      fcurrN->AddCross(oneSclP[i]);
716                      currentP->AddCross(oneSclN[tmplastNclIndex+k]);
717                      currentpkg->AddNSideCluster(oneSclN[lastNclIndex]);
718                      Nclpack++;
719                   }else break;
720               } else break;
721           }
722       
723           // end of package 
724           //if( IsCrossing )
725       } else { 
726          lastStripP  = currentP->GetLastDigitStripNo();
727          lastStripN  = currentN->GetLastDigitStripNo();
728          
729          if((lastStripN-fSFF) < lastStripP ) {
730             //Take new PCluster
731             if((Nclpack>0)&& (Pclpack>0)) {
732               new ((*fPackages)[fNPackages]) AliITSpackageSSD(fClusterP,fClusterN);
733               currentpkg = (AliITSpackageSSD*)((*fPackages)[fNPackages++]);
734             }
735               
736             //so we have to take another cluster on side N and check it
737             //we have to check until taken cluster's last strip will 
738             //be > last of this cluster(P)
739             //so it means that there is no corresponding cluster on side N 
740             //to this cluster (P)
741             //so we have to continue main loop with new "lastNclIndex"
742          
743             Nclpack=0;     
744             Pclpack=0;   
745             //We are not sure that next N cluter will cross with this P cluster
746             //There might one, or more, clusters that do not have any 
747             //corresponding cluster on the other side      
748             for(;;) {     
749                lastNclIndex++;
750                //Check if we are not out of array
751                if (lastNclIndex<fNClusterN) {
752                   currentN = GetNSideCluster(oneSclN[lastNclIndex]);
753                   if ( IsCrossing(currentP, currentN) ){
754                      Nclpack++;       
755                      Pclpack++;
756                      currentP->AddCross(oneSclN[lastNclIndex]);
757                      currentN->AddCross(oneSclP[i]);
758                      currentpkg->AddNSideCluster(oneSclN[lastNclIndex]);
759                      currentpkg->AddPSideCluster(oneSclP[i]);          
760                      //Check, if next N side clusters crosses with it too
761                      tmplastNclIndex=lastNclIndex;
762                      for(k=1;k<fSFF+1;k++) {
763                        //Check if we are not out of array
764                        if ((tmplastNclIndex+k)<fNClusterN) {
765                           fcurrN = GetNSideCluster(oneSclN[tmplastNclIndex+k]);
766                           if( IsCrossing(currentP, fcurrN) ) {
767                              Nclpack++;
768                              lastNclIndex++;
769                              currentP->AddCross(oneSclN[tmplastNclIndex+k]);
770                              fcurrN->AddCross(oneSclP[i]);
771                              currentpkg->AddNSideCluster(oneSclN[tmplastNclIndex+k]);
772                           }else break;
773                        } else break; //we are out of array
774                      } // end loop 
775                      break;
776                   } else {
777                      firstStripP = currentP->GetFirstDigitStripNo();
778                      firstStripN = currentN->GetFirstDigitStripNo();
779                      if((firstStripN+fSFB)>=firstStripP) break;
780                      else continue;
781                   }  
782                } else goto EndOfFunction;
783             } // end for(;;)
784          }  else  //EndOfPackage
785            {
786              continue;     
787            } //else EndOfPackage
788        
789       }//else IsCrossing
790             
791   }//main for
792
793  EndOfFunction:
794   if ((Nclpack<1) || (Pclpack<1)) fNPackages--;
795
796    delete oneSclP;
797    delete oneSclN;
798
799 }
800
801
802
803 //-----------------------------------------------
804 void AliITSClusterFinderSSD::PackagesToPoints()
805 {
806  
807  register Int_t i;
808  AliITSpackageSSD *currentpkg;
809  Int_t NumNcl;
810  Int_t NumPcl;
811  Int_t clusterIndex;
812  Bool_t clusterSide;
813
814  for(i=0;i<fNPackages;i++) {
815     //get pointer to next package
816     currentpkg = (AliITSpackageSSD*)((*fPackages)[i]); 
817     NumNcl = currentpkg->GetNumOfClustersN();
818     NumPcl = currentpkg->GetNumOfClustersP();
819     if(debug) cout<<"New Package\nNumPcl ="<<NumPcl<<" NumNcl ="<<NumNcl<<"\n";
820
821     while(((NumPcl>=2)&&(NumNcl>2))||((NumPcl>2)&&(NumNcl>=2))) {  
822         //This is case of "big" pacakage
823         //conditions see 3 lines above  
824         //if in the big package exists cluster with one cross 
825         //we can reconstruct this point without any geometrical ambiguities
826         if(currentpkg->GetClusterWithOneCross(clusterIndex, clusterSide) )
827         {
828             ResolveClusterWithOneCross(currentpkg, clusterIndex, clusterSide);
829         } else if (clusterIndex == -2) {
830             NumPcl = 0;
831             NumNcl = 0;
832             break;
833         } else {
834             if ( (NumNcl==NumPcl) && (NumPcl<10)) {
835                //if there is no cluster with one cross 
836                //we can resolve whole package at once 
837                //by finding best combination
838                //we can make combinations when NumNcl==NumPcl
839                if (ResolvePackageBestCombin(currentpkg)) { 
840                    //sometimes creating combinations fail, 
841                    //but it happens very rarely
842                    NumPcl = 0;
843                    NumNcl = 0;
844                    break;
845                } else ResolveOneBestMatchingPoint(currentpkg);
846             } else {
847                ResolveOneBestMatchingPoint(currentpkg);
848             }
849         }
850         NumNcl = currentpkg->GetNumOfClustersN();
851         NumPcl = currentpkg->GetNumOfClustersP();
852     } // end while 
853     if ((NumPcl==1)&&(NumNcl==1)) {
854       ResolveSimplePackage(currentpkg);
855        continue;
856     }
857     if (NumPcl==1) {
858       ResolvePackageWithOnePSideCluster(currentpkg);
859       continue;
860     }
861     if (NumNcl==1) {
862       ResolvePackageWithOneNSideCluster(currentpkg);
863       continue;
864     }
865     if ((NumPcl==2)&&(NumNcl==2)) {
866       ResolveTwoForTwoPackage(currentpkg);
867       continue; 
868     }
869     if ((NumPcl==0)&&(NumNcl>0)) { }
870     if ((NumNcl==0)&&(NumPcl>0)) { }
871
872   } // end loop over fNPackages
873
874
875 }
876
877
878 //----------------------------------------------------------
879
880 void AliITSClusterFinderSSD::
881 ResolveClusterWithOneCross(AliITSpackageSSD *currentpkg, Int_t clusterIndex, Bool_t clSide)
882 {
883
884    if (clSide == SIDEP) ResolvePClusterWithOneCross(currentpkg,clusterIndex);
885    else  ResolveNClusterWithOneCross(currentpkg,clusterIndex);
886
887 }
888
889
890 //---------------------------------------------------------
891 void AliITSClusterFinderSSD::
892 ResolvePClusterWithOneCross(AliITSpackageSSD *pkg, Int_t clusterIndex)
893 {
894
895 /*
896 There is cluster (side P) which crosses with only one cluster on the other 
897 side (N)
898
899 ie:
900     \    /  \/   /
901      \  /   /\  /
902       \/   /  \/
903       /\  /   /\
904      /  \/   /  \    .....
905     /   /\  /    \
906    /   /  \/      \
907   /   /   /\       \
908  /   /   /  \       \
909 */
910
911
912   AliITSclusterSSD * clusterP;
913   AliITSclusterSSD * clusterN;
914
915   Float_t posClusterP;           //Cluster P position in strip coordinates
916   Float_t posClusterN;           //Cluster N position in strip coordinates
917   
918   Float_t posErrorClusterP;
919   Float_t posErrorClusterN;
920   
921   Float_t sigClusterP;
922   Float_t sigClusterN;
923
924   Float_t sigErrorClusterP;
925   Float_t sigErrorClusterN;
926
927   Float_t ps;
928   Float_t ns;
929   
930   Float_t Chicomb;
931   Int_t clusterIdx;
932
933   if (debug) cout<<"ResolvePClusterWithOneCross\n";
934
935   clusterP=pkg->GetPSideCluster(clusterIndex);
936   posClusterP = GetClusterZ(clusterP);
937   posErrorClusterP = clusterP->GetPositionError();
938   sigClusterP = ps = clusterP->GetTotalSignal();
939   sigErrorClusterP = clusterP->GetTotalSignalError(); 
940   //carefully, it returns index in TClonesArray
941   
942   clusterIdx = clusterP->GetCross(0);
943   clusterN=this->GetNSideCluster(clusterIdx);
944   ns = clusterN->GetTotalSignal();
945   posClusterN = GetClusterZ(clusterN);
946   posErrorClusterN = clusterN->GetPositionError();
947   pkg->DelCluster(clusterIndex,SIDEP);
948   sigClusterN = ps/fPNsignalRatio;
949   // there is no sonse to check how signal ratio is far from perfect 
950   // matching line if the if below it is true
951   if (ns < sigClusterN) {
952       sigClusterN=ns;
953       if (debug) cout<<"n1 < p1/fPNsignalRatio";
954       if (debug) cout<<"Attempting to del cluster N "<<clusterIdx<<" ... \n";
955       pkg->DelClusterOI(clusterIdx,SIDEN);
956   } else {
957       //Let's see how signal ratio is far from perfect matching line
958       Chicomb = DistToPML(ps,ns);
959       if (debug) cout<<"Chic "<<Chicomb<<"\n";
960       if (Chicomb > falpha2) {
961          //it is near, so we can risk throwing this cluster away too
962          if (debug) cout<<"Attempting to del cluster N "<<clusterIdx<<"...\n"; 
963          pkg->DelClusterOI(clusterIdx,SIDEN);
964       } else {
965          clusterN->CutTotalSignal(sigClusterN);
966          if (debug) cout <<"Signal cut  |||||||||||||\n";
967       }
968   }
969   sigErrorClusterN= clusterN->GetTotalSignalError(); 
970   CreateNewRecPoint(posClusterP,posErrorClusterP,posClusterN,posErrorClusterN,
971                     sigClusterP+sigClusterN, sigErrorClusterN+sigErrorClusterP,
972                     clusterP, clusterN, 0.75);
973
974 }
975
976
977 //------------------------------------------------------
978 void AliITSClusterFinderSSD::
979 ResolveNClusterWithOneCross(AliITSpackageSSD *pkg, Int_t clusterIndex)
980 {
981
982   AliITSclusterSSD * clusterP;
983   AliITSclusterSSD * clusterN;
984
985   Float_t posClusterP;              //Cluster P position in strip coordinates
986   Float_t posClusterN;              //Cluster N position in strip coordinates
987
988   Float_t posErrorClusterP;
989   Float_t posErrorClusterN;
990
991   Float_t sigClusterP;
992   Float_t sigClusterN;
993
994   Float_t sigErrorClusterP;
995   Float_t sigErrorClusterN;
996
997   Float_t ps;
998   Float_t ns;
999   
1000   Float_t Chicomb;
1001   Int_t clusterIdx;
1002   
1003   if (debug) cout<<"ResolveNClusterWithOneCross\n"; 
1004   
1005   clusterN=pkg->GetNSideCluster(clusterIndex);
1006   posClusterN = GetClusterZ(clusterN);
1007   posErrorClusterN = clusterN->GetPositionError();
1008   sigClusterN = ns = clusterN->GetTotalSignal();
1009   sigErrorClusterN = clusterN->GetTotalSignalError(); 
1010   //carefully, it returns index in TClonesArray
1011   clusterIdx = clusterN->GetCross(0);
1012   clusterP=this->GetPSideCluster(clusterIdx);
1013   ps = clusterP->GetTotalSignal();
1014   posClusterP = GetClusterZ(clusterP);
1015   posErrorClusterP = clusterP->GetPositionError();
1016   pkg->DelCluster(clusterIndex,SIDEN);
1017   sigClusterP=ns*fPNsignalRatio;
1018   // there is no sonse to check how signal ratio is far from perfect 
1019   // matching line if the if below it is true
1020   if (ps < sigClusterP) {
1021       sigClusterP = ps;
1022       if (debug) cout<<"ps < ns*fPNsignalRatio";
1023       if (debug) cout<<"Attempting to del cluster P "<<clusterIdx<<" ... \n";
1024       pkg->DelClusterOI(clusterIdx,SIDEP);
1025   } else {
1026       //Let's see how signal ratio is far from perfect matching line
1027        Chicomb = DistToPML(ps,ns);
1028        if (debug) cout<<"Chic "<<Chicomb<<"\n";
1029        if (Chicomb > falpha2) {
1030           //it is near, so we can risk frowing this cluster away too
1031           if (debug) cout<<"Attempting to del cluster P "<<clusterIdx<<"...\n";
1032           pkg->DelClusterOI(clusterIdx,SIDEP);
1033        } else {
1034           clusterN->CutTotalSignal(sigClusterP);
1035           if (debug) cout <<"Signal cut  ------------\n";
1036        }
1037   }
1038   sigErrorClusterP= clusterP->GetTotalSignalError(); 
1039   CreateNewRecPoint( posClusterP,posErrorClusterP,posClusterN,posErrorClusterN,
1040                      sigClusterP+sigClusterN,sigErrorClusterN+sigErrorClusterP,
1041                      clusterP, clusterN, 0.75);
1042
1043
1044 }
1045
1046
1047
1048 //-------------------------------------------------------
1049
1050 Bool_t AliITSClusterFinderSSD::
1051 ResolvePackageBestCombin(AliITSpackageSSD *pkg)
1052 {
1053   if (debug) cout<<"NumNcl==NumPcl ("<<pkg->GetNumOfClustersN()
1054       <<"=="<<pkg->GetNumOfClustersP()<<"); Generating combinations ... \n";
1055
1056
1057   AliITSclusterSSD * clusterP;
1058   AliITSclusterSSD * clusterN;
1059
1060   Int_t Ncombin; //Number of combinations
1061   Int_t itera;   //iterator
1062   Int_t sizet=1;   //size of array to allocate 
1063  
1064   Int_t NP = pkg->GetNumOfClustersP();
1065   for(itera =2; itera <= NP ;itera ++) {
1066      sizet=sizet*itera;
1067      if (sizet > 10000) {
1068        sizet=10000;
1069        break;
1070      }
1071   }
1072
1073   Int_t **combin = new Int_t*[sizet]; //2D array to keep combinations in
1074
1075   for(itera =0; itera <sizet;itera++) {
1076    combin[itera] = new Int_t[NP+1];
1077   }
1078      
1079   pkg->GetAllCombinations(combin,Ncombin,sizet);
1080  
1081   if (Ncombin==0) { 
1082     if (debug) cout<<"No combin Error";
1083     return kFALSE;
1084   }
1085 //calculate which combination fits the best to perfect matching line 
1086   Int_t bc = GetBestComb(combin,Ncombin,NP,pkg); 
1087   if (debug) cout<<"  bc "<<bc <<"\n";
1088               
1089   for(itera =0; itera < NP; itera ++) {
1090       clusterP = pkg->GetPSideCluster(itera);
1091       //carefully here
1092       //becase AliITSclusterSSD::GetCross returns index in 
1093       //AliITSmoduleSSD.fClusterP, which is put to "combin"        
1094       clusterN = GetNSideCluster(combin[bc][itera]);
1095       CreateNewRecPoint(clusterP, clusterN, 0.75);
1096   } 
1097
1098   for(itera =0; itera <sizet;itera++) {
1099      delete [](combin[itera]);
1100   }
1101   delete [] combin;  
1102   return kTRUE;
1103
1104 }
1105
1106
1107 //----------------------------------------------------
1108 void AliITSClusterFinderSSD::
1109 ResolveOneBestMatchingPoint(AliITSpackageSSD *pkg)
1110 {
1111
1112  Int_t ni, pi;
1113
1114  Int_t prvP, prvN;
1115  Int_t nextP, nextN=0;
1116
1117   
1118  AliITSclusterSSD * clusterP;
1119  AliITSclusterSSD * clusterN;
1120
1121  Bool_t split = kFALSE;
1122
1123  if (debug) cout<<"ResolveOneBestMatchingPoint \n";
1124
1125  GetBestMatchingPoint(pi, ni, pkg);
1126          
1127  clusterP = GetPSideCluster(pi);
1128  clusterN = GetNSideCluster(ni);
1129               
1130  CreateNewRecPoint(clusterP, clusterN, 0.75);
1131  
1132  if ((nextP=pkg->GetNextPIdx(pi))!=-1)
1133    if ((nextN=pkg->GetNextNIdx(ni))!=-1)
1134      if ((prvP=pkg->GetPrvPIdx(pi))!=-1)
1135        if ((prvN=pkg->GetPrvNIdx(ni))!=-1)
1136           if( !(GetPSideCluster(prvP)->IsCrossingWith(nextN)))
1137             if( !(GetPSideCluster(nextP)->IsCrossingWith(prvN)))
1138               {
1139                  split = kTRUE;
1140               }   
1141  
1142  
1143  pkg->DelClusterOI(pi, SIDEP);
1144  pkg->DelClusterOI(ni, SIDEN);
1145  
1146  if (split) {
1147    if (debug) cout<<"spltting package ...\n";
1148    new ((*fPackages)[fNPackages]) AliITSpackageSSD(fClusterP,fClusterN);
1149    pkg->
1150      SplitPackage(nextP,nextN,(AliITSpackageSSD*)((*fPackages)[fNPackages++]));
1151  }
1152
1153 }
1154
1155
1156 //--------------------------------------------------
1157 void  AliITSClusterFinderSSD::ResolveSimplePackage(AliITSpackageSSD *pkg)
1158 {
1159
1160   AliITSclusterSSD * clusterP;
1161   AliITSclusterSSD * clusterN;
1162
1163   clusterP         = pkg->GetPSideCluster(0);
1164         
1165   clusterN         = pkg->GetNSideCluster(0);
1166
1167   CreateNewRecPoint(clusterP, clusterN, 1.0);
1168
1169
1170
1171
1172
1173 //--------------------------------------------------
1174 void AliITSClusterFinderSSD:: ResolvePackageWithOnePSideCluster(AliITSpackageSSD *pkg) {
1175
1176
1177 /*
1178  \   \   \  /
1179   \   \   \/
1180    \   \  /\
1181     \   \/  \
1182      \  /\   \
1183       \/  \   \
1184       /\   \   \
1185      /  \   \   \  
1186
1187 */
1188
1189
1190 /************************/
1191 /**************************/
1192 /*XP SP itg jest tylko jeden nie musi byc tablica i w peetli nie trzeba po niej iterowcac*/
1193 /***************************/
1194
1195
1196  Int_t k;
1197  AliITSclusterSSD * clusterP;
1198  AliITSclusterSSD * clusterN;
1199  
1200  Int_t NN=pkg->GetNumOfClustersN();
1201  Float_t sumsig = 0;
1202  
1203  Float_t   XP;
1204  Float_t   XPerr;
1205  
1206  Float_t * XN = new Float_t[NN];
1207  Float_t * XNerr = new Float_t[NN];
1208  
1209  Float_t * SP = new Float_t[NN];
1210  Float_t * SN = new Float_t[NN];
1211  
1212  Float_t * SPerr = new Float_t[NN];
1213  Float_t * SNerr = new Float_t[NN];
1214  
1215  Float_t p1;
1216  Float_t p1err;
1217  
1218  clusterP = pkg->GetPSideCluster(0);
1219  p1       = clusterP->GetTotalSignal();
1220    
1221  XP    = GetClusterZ(clusterP);
1222  XPerr = clusterP->GetPositionError();
1223  p1err = clusterP->GetTotalSignalError();
1224    
1225  for(k=0;k<NN;k++) {
1226     SN[k]    = pkg->GetNSideCluster(k)->GetTotalSignal();
1227     SNerr[k] = pkg->GetNSideCluster(k)->GetTotalSignalError();
1228     sumsig   += SN[k];
1229  }
1230  for(k=0;k<NN;k++) {
1231     clusterN = pkg->GetNSideCluster(k);
1232     SP[k]= p1*SN[k]/sumsig;
1233     SPerr[k] = p1err*SN[k]/sumsig;
1234     XN[k]=GetClusterZ(clusterN);
1235     XNerr[k]= clusterN->GetPositionError();
1236     CreateNewRecPoint(XP,XPerr, XN[k],XNerr[k], SP[k]+SN[k], SPerr[k]+SNerr[k], clusterP, clusterN, 1.0);
1237
1238  }
1239
1240 }
1241
1242
1243 //---------------------------------------------------------
1244 void AliITSClusterFinderSSD::ResolvePackageWithOneNSideCluster(AliITSpackageSSD *pkg) {
1245
1246
1247 /*
1248     \    /   /   /
1249      \  /   /   /
1250       \/   /   /
1251       /\  /   /
1252      /  \/   /
1253     /   /\  /  
1254    /   /  \/
1255   /   /   /\
1256  /   /   /  \  
1257 */
1258
1259
1260  Int_t k;
1261  AliITSclusterSSD * clusterP;
1262  AliITSclusterSSD * clusterN;
1263  
1264  Int_t NP=pkg->GetNumOfClustersP();
1265  Float_t sumsig = 0;
1266  
1267  Float_t * XP = new Float_t[NP];
1268  Float_t * XPerr = new Float_t[NP];
1269  
1270  Float_t   XN;
1271  Float_t   XNerr;
1272  
1273  Float_t * SP = new Float_t[NP];
1274  Float_t * SN = new Float_t[NP];
1275  
1276  Float_t * SPerr = new Float_t[NP];
1277  Float_t * SNerr = new Float_t[NP];
1278  
1279  Float_t n1;
1280  Float_t n1err;
1281
1282  clusterN = pkg->GetNSideCluster(0);
1283  n1       = clusterN->GetTotalSignal();
1284    
1285  XN = GetClusterZ(clusterN);
1286  XNerr = clusterN->GetPositionError();
1287  
1288  n1err=clusterN->GetTotalSignalError();
1289    
1290  for(k=0;k<NP;k++) {
1291     SP[k] = pkg->GetPSideCluster(k)->GetTotalSignal();
1292     sumsig += SP[k];
1293     SPerr[k] = pkg->GetPSideCluster(k)->GetTotalSignalError();
1294  }
1295
1296  for(k=0;k<NP;k++) {
1297     clusterP = pkg->GetPSideCluster(k);
1298     SN[k]= n1*SP[k]/sumsig;
1299     XP[k]=GetClusterZ(clusterP);
1300     XPerr[k]= clusterP->GetPositionError();
1301     SNerr[k] = n1err*SP[k]/sumsig;
1302     CreateNewRecPoint(XP[k],XPerr[k], XN,XNerr, SP[k]+SN[k], SPerr[k]+SNerr[k],clusterP, clusterN, 1.0);
1303   }
1304   
1305                
1306 }
1307
1308 /**********************************************************************/
1309 /*********      2  X   2  *********************************************/
1310 /**********************************************************************/
1311
1312 void AliITSClusterFinderSSD::
1313 ResolveTwoForTwoPackage(AliITSpackageSSD *pkg)
1314 {
1315
1316   AliITSclusterSSD *clusterP1 = pkg->GetPSideCluster(0);
1317   AliITSclusterSSD *clusterP2 = pkg->GetPSideCluster(1);
1318   AliITSclusterSSD *clusterN1 = pkg->GetNSideCluster(0);
1319   AliITSclusterSSD *clusterN2 = pkg->GetNSideCluster(1);
1320
1321   AliITSclusterSSD *tmp;
1322
1323   Float_t p1sig;
1324   Float_t p2sig;
1325   Float_t n1sig;
1326   Float_t n2sig;
1327
1328   Float_t p1sigErr;
1329   Float_t p2sigErr;
1330   Float_t n1sigErr;
1331   Float_t n2sigErr;
1332
1333   Float_t ZposP1;
1334   Float_t ZposP2;
1335   Float_t ZposN1;
1336   Float_t ZposN2;
1337
1338   Float_t ZposErrP1;
1339   Float_t ZposErrP2;
1340   Float_t ZposErrN1;
1341   Float_t ZposErrN2;
1342
1343   Float_t D12;
1344
1345   Float_t Chicomb1;
1346   Float_t Chicomb2;
1347
1348   if(clusterP1->GetDigitStripNo(0) > clusterP2->GetDigitStripNo(0)) {
1349      if (debug) cout<<"P strips flip\n";
1350      tmp       = clusterP1;
1351      clusterP1 = clusterP2;
1352      clusterP2 = tmp;
1353   }
1354   if(clusterN1->GetDigitStripNo(0) > clusterN2->GetDigitStripNo(0)) {
1355      if (debug) cout<<"N strips flip\n";
1356      tmp       = clusterN1;
1357      clusterN1 = clusterN2;
1358      clusterN2 = tmp;
1359   }
1360   p1sig = clusterP1->GetTotalSignal();
1361   p2sig = clusterP2->GetTotalSignal();
1362   n1sig = clusterN1->GetTotalSignal();
1363   n2sig = clusterN2->GetTotalSignal();
1364
1365
1366   p1sigErr = clusterP1->GetTotalSignalError();
1367   n1sigErr = clusterN1->GetTotalSignalError();
1368   p2sigErr = clusterP2->GetTotalSignalError();
1369   n2sigErr = clusterN2->GetTotalSignalError();
1370           
1371   ZposP1 = clusterP1->GetPosition();
1372   ZposN1 = clusterN1->GetPosition();
1373   ZposP2 = clusterP2->GetPosition();
1374   ZposN2 = clusterN2->GetPosition();
1375
1376   ZposErrP1 = clusterP1->GetPositionError();
1377   ZposErrN1 = clusterN1->GetPositionError();
1378   ZposErrP2 = clusterP2->GetPositionError();
1379   ZposErrN2 = clusterN2->GetPositionError();
1380
1381   //in this may be two types:
1382   //    1.When each cluster crosses with 2 clusters on the other side
1383   //            gives 2 ghosts and 2 real points
1384   //
1385   //    2.When two clusters (one an each side) crosses with only one on 
1386   //           the other side and two crosses (one on the each side) with 
1387   //           two on the other gives 2 or 3 points,
1388    
1389   if (debug) cout<<"2 for 2 ambiguity ...";
1390   /***************************/
1391   /*First sort of combination*/
1392   /***************************/
1393   
1394   if((clusterP1->GetCrossNo()==2) && (clusterN1->GetCrossNo()==2)) {  
1395      if (debug) cout<<"All clusters has 2 crosses\n";
1396      D12 = TMath::Sqrt((Float_t)((p1sig -p2sig)*(p1sig -p2sig) + (n1sig -n2sig)*(n1sig -n2sig)));
1397             
1398      if(debug) cout<<"Distance between points in P(N) plane D12 = "<<D12<<"\n";
1399             
1400      /*********************************************/
1401      /*resolving ambiguities:                     */
1402      /*we count Chicomb's and we take combination */
1403      /*giving lower Chicomb as a real points      */
1404      /*Keep only better combinantion              */
1405      /*********************************************/
1406             
1407      if (D12 > (falpha3*17.8768)) {
1408        if (debug) cout<<"decided to take only one pair \n";
1409        Chicomb1 = DistToPML(p1sig,n1sig) + DistToPML(p2sig,n2sig);
1410        Chicomb2 = DistToPML(p2sig,n1sig) + DistToPML(p1sig,n2sig);
1411        if (debug) {
1412          cout<<" 00 11 combination : "<<Chicomb1<<" 01 10 combination : "<<Chicomb2<<" \n";
1413          cout<<"p1 = "<<p1sig<<"  n1 = "<<n1sig<<"  p2 = "<<p2sig<<"  n2 = "<<n2sig<<"\n"; 
1414        }
1415        if (Chicomb1 < Chicomb2) {
1416           if (debug) cout<<"00 11";
1417           CreateNewRecPoint(ZposP1,ZposErrP1, ZposN1,ZposErrN1, p1sig+n1sig, p1sigErr+n1sigErr, clusterP1, clusterN1, 0.75);   
1418           CreateNewRecPoint(ZposP2,ZposErrP2, ZposN2,ZposErrN2, p2sig+n2sig, p2sigErr+n2sigErr, clusterP2, clusterN2, 0.75); 
1419
1420   
1421           //second cominantion
1422        } else {
1423           if (debug) cout<<"01 10";
1424           CreateNewRecPoint(ZposP1,0, ZposN2,0, p1sig+n2sig, p1sigErr+n2sigErr, clusterP1, clusterN2, 0.75);   
1425           CreateNewRecPoint(ZposP2,0, ZposN1,0, p2sig+n1sig, p2sigErr+n1sigErr, clusterP2, clusterN1, 0.75);    
1426        } //end second combinantion
1427        //if (D12 > falpha3*17.8768)
1428        //keep all combinations
1429      } else {
1430         if (debug) cout<<"We decide to take all points\n";
1431         CreateNewRecPoint(ZposP1,ZposErrP1, ZposN1,ZposErrN1, p1sig+n1sig, p1sigErr+n1sigErr, clusterP1, clusterN1, 0.5);   
1432         CreateNewRecPoint(ZposP2,ZposErrP2, ZposN2,ZposErrN2, p2sig+n2sig, p2sigErr+n2sigErr, clusterP2, clusterN2, 0.5);   
1433         CreateNewRecPoint(ZposP1,ZposErrP1, ZposN2,ZposErrN2, p1sig+n2sig, p1sigErr+n2sigErr, clusterP1, clusterN2, 0.5);   
1434         CreateNewRecPoint(ZposP2,ZposErrP2, ZposN1,ZposErrN1, p2sig+n1sig, p2sigErr+n1sigErr, clusterP2, clusterN1, 0.5);    
1435      }
1436   }
1437   // ad.2 Second type of combination
1438   else {
1439     Chicomb1 = DistToPML(p1sig,n1sig) + DistToPML(p2sig,n2sig);
1440     if (debug) cout<<"\nhere can be reconstructed 3 points: chicomb = "<<Chicomb1<<"\n"; 
1441   
1442     if (Chicomb1<falpha1) {
1443        if (debug) cout<<"\nWe decided to take 3rd point"; 
1444        if (clusterP1->GetCrossNo()==1) {
1445           if (debug) cout<<"...  P1 has one cross\n"; 
1446           n1sig = p1sig/fPNsignalRatio;
1447           p2sig = n2sig*fPNsignalRatio;
1448               
1449           clusterN1->CutTotalSignal(n1sig);
1450           clusterP2->CutTotalSignal(p2sig);
1451        
1452           CreateNewRecPoint(ZposP2,ZposErrP2, ZposN1,ZposErrN1, p2sig+n1sig, p2sigErr+n1sigErr, clusterP2, clusterN1, 0.5);    
1453
1454           n1sig = clusterN1->GetTotalSignal();
1455           p2sig = clusterP2->GetTotalSignal();
1456
1457        } else {
1458           if (debug) cout<<"...  N1 has one cross\n";
1459               
1460           n2sig=p2sig/fPNsignalRatio;
1461           p1sig=n1sig*fPNsignalRatio;
1462
1463           clusterN2->CutTotalSignal(n2sig);
1464           clusterP1->CutTotalSignal(p1sig);
1465
1466           CreateNewRecPoint(ZposP1,ZposErrP1, ZposN2,ZposErrN2, p1sig+n2sig, p1sigErr+n2sigErr, clusterP1, clusterN2, 0.5);   
1467           
1468           n2sig=clusterN2->GetTotalSignal();
1469           p1sig=clusterP1->GetTotalSignal();
1470        }
1471     } else {
1472        if (debug) cout<<"\nWe decided  NOT to take 3rd point\n"; 
1473     }
1474
1475     CreateNewRecPoint(ZposP1,ZposErrP1, ZposN1,ZposErrN1, p1sig+n1sig, p1sigErr+n1sigErr, clusterP1, clusterN1, 1.0);   
1476     CreateNewRecPoint(ZposP2,ZposErrP2, ZposN2,ZposErrN2, p2sig+n2sig, p2sigErr+n2sigErr, clusterP2, clusterN2, 1.0);   
1477   
1478   }
1479
1480 }
1481
1482
1483
1484 //------------------------------------------------------
1485 Bool_t AliITSClusterFinderSSD::    
1486 CreateNewRecPoint(Float_t P, Float_t dP, Float_t N, Float_t dN,
1487                   Float_t Sig,Float_t dSig, 
1488                   AliITSclusterSSD *clusterP, AliITSclusterSSD *clusterN,
1489                   Stat_t prob)
1490 {
1491
1492   const Float_t kdEdXtoQ = 2.778e+8;
1493   const Float_t kconv = 1.0e-4; 
1494   const Float_t kRMSx = 20.0*kconv; // microns->cm ITS TDR Table 1.3
1495   const Float_t kRMSz = 830.0*kconv; // microns->cm ITS TDR Table 1.3
1496
1497
1498   Float_t p=P;
1499   Float_t n=N;
1500   if (GetCrossing(P,N)) {
1501      GetCrossingError(dP,dN);
1502      AliITSRawClusterSSD cnew;
1503      Int_t nstripsP=clusterP->GetNumOfDigits();
1504      Int_t nstripsN=clusterN->GetNumOfDigits();
1505      cnew.fMultiplicity=nstripsP;
1506      cnew.fMultiplicityN=nstripsN;
1507      /*
1508      if (nstripsP > 100) {
1509         printf("multiplicity > 100 - increase the dimension of the arrays %d\n",nstripsP);
1510         nstripsP=100;
1511      }
1512      Int_t i;
1513      for(i=0;i<nstripsP;i++) {
1514        // check if 'clusterP->GetDigitStripNo(i)' returns the digit index
1515        cnew.fIndexMap[i] = clusterP->GetDigitStripNo(i); 
1516      } 
1517      if (nstripsN > 100) {
1518         printf("multiplicity > 100 - increase the dimension of the arrays %d\n",nstripsN);
1519         nstripsN=100;
1520      }
1521      for(i=0;i<nstripsN;i++) {
1522        // check if 'clusterN->GetDigitStripNo(i)' returns the digit index
1523        cnew.fIndexMapN[i] = clusterN->GetDigitStripNo(i); 
1524      }
1525      */
1526      cnew.fQErr=dSig;
1527      //cnew.fProbability=(float)prob; 
1528      fITS->AddCluster(2,&cnew);
1529      // add the rec point info
1530      AliITSRecPoint rnew;
1531      rnew.SetX(P*kconv);
1532      rnew.SetZ(N*kconv);
1533      rnew.SetQ(Sig);
1534      rnew.SetdEdX(Sig/kdEdXtoQ);
1535      rnew.SetSigmaX2( kRMSx* kRMSx); 
1536      rnew.SetSigmaZ2( kRMSz* kRMSz);
1537      rnew.SetProbability((float)prob);
1538      fITS->AddRecPoint(rnew);
1539      /*
1540      // it was
1541      fPointsM->AddLast( (TObject*) 
1542                    new AliITSRecPointSSD(P,dP,0,0,N,dN,Sig,dSig,fLayer,fLad,fDet,clusterP,clusterN,prob) );
1543      */
1544     
1545      if(debug) cout<<"\n"<<": ImpactPoint Created: X = "<<P<<" Z = "<<N<<";   P = "<<p<<"  N = "<<n<<"\n";
1546      return kTRUE;
1547   } else { 
1548      if (debug) {
1549        cout<<"\n"<<": ATENTION : stupid ImpactPoit Point X = "<<P<<" Z = "<<N<<";  P = "<<p<<"  N = "<<n<<"\n";
1550      }
1551      return kFALSE;
1552   }
1553 }
1554
1555
1556
1557 /**********************************************************************/
1558 Bool_t  AliITSClusterFinderSSD::
1559 CreateNewRecPoint(AliITSclusterSSD *clusterP, AliITSclusterSSD *clusterN, Stat_t prob)
1560 {
1561   Float_t posClusterP;  //Cluster P position in strip coordinates
1562   Float_t posClusterN;  //Cluster N position in strip coordinates
1563  
1564   Float_t posErrorClusterP;
1565   Float_t posErrorClusterN;
1566
1567   Float_t sigClusterP;
1568   Float_t sigClusterN;
1569
1570   Float_t sigErrorClusterP;
1571   Float_t sigErrorClusterN;
1572
1573   posClusterP      = clusterP->GetPosition();
1574   posErrorClusterP = clusterP->GetPositionError();
1575   sigClusterP      = clusterP->GetTotalSignal();
1576   sigErrorClusterP = clusterP->GetTotalSignalError();
1577
1578   posClusterN      = clusterN->GetPosition();
1579   posErrorClusterN = clusterN->GetPositionError();
1580   sigClusterN      = clusterN->GetTotalSignal();
1581   sigErrorClusterN = clusterN->GetTotalSignalError();
1582
1583   return CreateNewRecPoint( posClusterP,posErrorClusterP,  posClusterN,posErrorClusterN, 
1584                      sigClusterP+sigClusterN, sigErrorClusterN+sigErrorClusterP,
1585                      clusterP, clusterN, prob);
1586
1587 }
1588
1589 //--------------------------------------------------
1590 Bool_t AliITSClusterFinderSSD::IsCrossing(AliITSclusterSSD* p, AliITSclusterSSD* n)
1591 {
1592   Float_t x = p->GetPosition();
1593   Float_t y = n->GetPosition();
1594   return GetCrossing(x,y);
1595 }
1596
1597
1598 //----------------------------------------------
1599 Bool_t AliITSClusterFinderSSD::Strip2Local( Float_t stripP, Float_t stripN, Float_t &Z,Float_t &X)
1600 {
1601
1602
1603 //Piotr Krzysztof Skowronski
1604 //Warsaw University of Technology
1605 //skowron@if.pw.edu.pl
1606
1607 /*
1608  Z = (stripN-stripP)*fFactorOne;  
1609  X = (stripN + stripP - fStripZeroOffset)*fFactorTwo;
1610 */
1611
1612  Float_t P =stripP;
1613  Float_t N =stripN;
1614
1615  GetCrossing(P,N);
1616  X=P;
1617  Z=N;
1618  if (debug) cout<<"P="<<stripP<<" N="<<stripN<<"   X = "<<X<<" Z = "<<Z<<"\n";
1619  if ((Z<2.1)&&(Z>-2.1)&&(X<3.65)&&(X>-3.65)) return true; 
1620     else return false;  
1621
1622 }
1623
1624
1625 //-----------------------------------------------------------
1626 Float_t  AliITSClusterFinderSSD::GetClusterZ(AliITSclusterSSD* clust)
1627 {
1628
1629   return clust->GetPosition();
1630
1631 }
1632
1633 //-------------------------------------------------------------
1634 Int_t AliITSClusterFinderSSD::GetBestComb
1635 (Int_t** comb,Int_t Ncomb, Int_t Ncl, AliITSpackageSSD * pkg)
1636 {
1637
1638
1639 //Piotr Krzysztof Skowronski
1640 //Warsaw University of Technology
1641 //skowron@if.pw.edu.pl
1642
1643 //returns index of best combination in "comb"
1644 //comb : sets of combinations
1645 //       in given combination on place "n" is index of 
1646 //       Nside cluster coresponding to "n" P side cluster
1647 //
1648 //Ncomb : number of combinations == number of rows in "comb"
1649 //
1650 //Ncl   : number of clusters in each row  == number of columns in "comb"
1651 //
1652 //pkg   : package 
1653
1654
1655   Float_t currbval=-1;  //current best value, 
1656                         //starting value is set to number negative, 
1657                         //to be changed in first loop turn automatically
1658       
1659   Int_t  curridx=0;     //index of combination giving currently the best chi^2 
1660   Float_t chi;
1661   Float_t ps, ns;       //signal of P cluster and N cluster
1662
1663   Int_t i,j;
1664   for(i=0;i<Ncomb;i++)
1665     {
1666       chi=0;
1667       
1668       for(j=0;j<Ncl;j++)
1669       {
1670         ps = pkg->GetPSideCluster(j)->GetTotalSignal();  //carrefully here, different functions
1671         ns = GetNSideCluster(comb[i][j])->GetTotalSignal(); 
1672         chi+=DistToPML(ps,ns);
1673       }
1674       
1675       if (currbval==-1) currbval = chi;
1676       if (chi<currbval)
1677         {
1678           curridx = i;
1679           currbval  =chi;
1680         }
1681     }
1682
1683
1684   return curridx;
1685
1686 }
1687
1688 //--------------------------------------------------
1689 void AliITSClusterFinderSSD::GetBestMatchingPoint
1690 (Int_t & ip, Int_t & in, AliITSpackageSSD* pkg )
1691 {
1692
1693
1694 //Piotr Krzysztof Skowronski
1695 //Warsaw University of Technology
1696 //skowron@if.pw.edu.pl
1697  
1698   if (debug) pkg->PrintClusters();
1699  
1700   Float_t ps, ns;  //signals on side p & n
1701   
1702   Int_t nc;       // number of crosses in given p side cluster
1703   Int_t n,p;
1704   AliITSclusterSSD *curPcl, *curNcl;  //current p side cluster
1705   Float_t  bestchi, chi;
1706   ip=p=pkg->GetPSideClusterIdx(0);
1707   in=n=pkg->GetNSideClusterIdx(0);
1708         
1709   bestchi=DistToPML(  pkg->GetPSideCluster(0)->GetTotalSignal(),
1710                   pkg->GetNSideCluster(0)->GetTotalSignal()  );
1711   Int_t i,j;
1712   for(i = 0; i< pkg->GetNumOfClustersP(); i++)
1713     {
1714       
1715       p =  pkg->GetPSideClusterIdx(i);   
1716       curPcl= GetPSideCluster(p);
1717       nc=curPcl->GetCrossNo();
1718       ps=curPcl->GetTotalSignal();
1719       
1720       for(j = 0; j< nc; j++)
1721         {
1722          n=curPcl->GetCross(j);
1723          curNcl= GetNSideCluster(n);
1724          ns=curNcl->GetTotalSignal();
1725          chi = DistToPML(ps, ns);
1726     
1727          if (chi <= bestchi)
1728           {
1729             bestchi = chi;
1730             in = n;
1731             ip = p;
1732           }
1733         }
1734     }
1735   
1736 }
1737
1738
1739 //------------------------------------------------------
1740 void  AliITSClusterFinderSSD::CalcStepFactor(Float_t Psteo, Float_t Nsteo )
1741 {
1742
1743
1744 //Piotr Krzysztof Skowronski
1745 //Warsaw University of Technology
1746 //skowron@if.pw.edu.pl
1747
1748
1749   // 95 is the pitch, 4000 - dimension along z ?
1750   /*
1751   fSFF = ( (Int_t)  (Psteo*40000/95 ) );// +1;
1752   fSFB = ( (Int_t)  (Nsteo*40000/95 ) );// +1;
1753   */
1754
1755
1756   Float_t dz=fSegmentation->Dz();
1757
1758   fSFF = ( (Int_t)  (Psteo*dz/fPitch ) );// +1;
1759   fSFB = ( (Int_t)  (Nsteo*dz/fPitch ) );// +1;
1760
1761 }
1762
1763
1764 //-----------------------------------------------------------
1765 AliITSclusterSSD* AliITSClusterFinderSSD::GetPSideCluster(Int_t idx)
1766 {
1767
1768
1769 //Piotr Krzysztof Skowronski
1770 //Warsaw University of Technology
1771 //skowron@if.pw.edu.pl
1772
1773
1774
1775   if((idx<0)||(idx>=fNClusterP))
1776     {
1777       printf("AliITSClusterFinderSSD::GetPSideCluster  : index out of range\n");
1778       return 0;
1779     }
1780   else
1781     {
1782       return (AliITSclusterSSD*)((*fClusterP)[idx]);
1783     }
1784 }
1785
1786 //-------------------------------------------------------
1787 AliITSclusterSSD* AliITSClusterFinderSSD::GetNSideCluster(Int_t idx)
1788 {
1789
1790
1791 //Piotr Krzysztof Skowronski
1792 //Warsaw University of Technology
1793 //skowron@if.pw.edu.pl
1794
1795   if((idx<0)||(idx>=fNClusterN))
1796     {
1797       printf("AliITSClusterFinderSSD::GetNSideCluster  : index out of range\n");
1798       return 0;
1799     }
1800   else
1801     {
1802       return (AliITSclusterSSD*)((*fClusterN)[idx]);
1803     }
1804 }
1805
1806 //--------------------------------------------------------
1807 AliITSclusterSSD* AliITSClusterFinderSSD::GetCluster(Int_t idx, Bool_t side)
1808 {
1809
1810
1811 //Piotr Krzysztof Skowronski
1812 //Warsaw University of Technology
1813 //skowron@if.pw.edu.pl
1814
1815   return (side) ? GetPSideCluster(idx) : GetNSideCluster(idx);
1816 }
1817
1818
1819 //--------------------------------------------------------
1820 void AliITSClusterFinderSSD::ConsumeClusters()
1821 {
1822
1823  Int_t i;
1824  for(i=0;i<fNPackages;i++)
1825   {
1826     ((AliITSpackageSSD*)((*fPackages)[i]))->ConsumeClusters();
1827   }
1828
1829 }
1830
1831
1832 //--------------------------------------------------------
1833 void AliITSClusterFinderSSD::ReconstructNotConsumedClusters()
1834 {
1835   Int_t i;
1836   AliITSclusterSSD *cluster;
1837   Float_t pos;
1838   Float_t sig;
1839   Float_t sigerr;
1840   Float_t x1,x2,z1,z2;
1841
1842   Float_t Dx = fSegmentation->Dx();
1843   Float_t Dz = fSegmentation->Dz();
1844   
1845   const Float_t kdEdXtoQ = 2.778e+8; // GeV -> number of e-hole pairs
1846   const Float_t kconv = 1.0e-4; 
1847   const Float_t kRMSx = 20.0*kconv; // microns->cm ITS TDR Table 1.3
1848   const Float_t kRMSz = 830.0*kconv; // microns->cm ITS TDR Table 1.3
1849   
1850   for(i=0;i<fNClusterP;i++)
1851     {
1852       //printf("P side cluster: cluster number %d\n",i);
1853       cluster = GetPSideCluster(i);
1854       if (!cluster->IsConsumed())
1855         {
1856           if( (pos = cluster->GetPosition()) < fSFB)
1857             {
1858               sig = cluster->GetTotalSignal();
1859               
1860               sig += sig/fPNsignalRatio;
1861               
1862               sigerr = cluster->GetTotalSignalError();
1863               x1 = -Dx/2 + pos *fPitch;
1864               z1 = -Dz/2;
1865
1866               x2 = pos;
1867               z2 = 1;
1868               GetCrossing (x2,z2);
1869               AliITSRawClusterSSD cnew;
1870               Int_t nstripsP=cluster->GetNumOfDigits();
1871               cnew.fMultiplicity=nstripsP;
1872               cnew.fMultiplicityN=0;
1873               /*
1874               if (nstripsP > 100) {
1875                 printf("multiplicity > 100 - increase the dimension of the arrays %d\n",nstripsP);
1876                 nstripsP=100;
1877               }
1878               Int_t k;
1879               for(k=0;k<nstripsP;k++) {
1880                 // check if 'clusterP->GetDigitStripNo(i)' returns 
1881                 // the digit index and not smth else
1882                 cnew.fIndexMap[k] = cluster->GetDigitStripNo(k); 
1883               } 
1884               */
1885               cnew.fQErr=sigerr;
1886               //cnew.fProbability=0.75; 
1887               fITS->AddCluster(2,&cnew);
1888               // add the rec point info
1889               AliITSRecPoint rnew;
1890               rnew.SetX(kconv*(x1+x2)/2);
1891               rnew.SetZ(kconv*(z1+z2)/2);
1892               rnew.SetQ(sig);
1893               rnew.SetdEdX(sig/kdEdXtoQ);
1894               rnew.SetSigmaX2( kRMSx* kRMSx);
1895               rnew.SetSigmaZ2( kRMSz* kRMSz);
1896               rnew.SetProbability(0.75);
1897               fITS->AddRecPoint(rnew);
1898               /*
1899               fPointsM->AddLast( (TObject*) 
1900                    new AliITSRecPointSSD((x1+x2)/2 ,0,0,0,(z1+z2)/2,0,sig,sigerr,fLayer,fLad,fDet,cluster, 0.75) );
1901               */
1902
1903               if(debug) cout<<"\n"<<": SINGLE SIDE ImpactPoint Created: X = "
1904                   <<(x1+x2)/2<<" Z = "<<(z1+z2)/2<<";   Pos = "<<pos;
1905             }
1906         }
1907     }
1908
1909   for(i=0;i<fNClusterN;i++)
1910     {
1911       // printf("N side cluster: cluster number %d\n",i);
1912       cluster = GetNSideCluster(i);
1913       if (!cluster->IsConsumed())
1914         {
1915           if( (pos = cluster->GetPosition()) < fSFF)
1916             {
1917               sig = cluster->GetTotalSignal();
1918               
1919               sig += sig*fPNsignalRatio;
1920               
1921               sigerr = cluster->GetTotalSignalError();
1922               
1923               x1 = -Dx/2 + pos *fPitch;
1924               z1 = Dz/2;
1925
1926               x2 = 1;
1927               z2 = pos;
1928               
1929               GetCrossing (x2,z2);
1930               AliITSRawClusterSSD cnew;
1931               Int_t nstripsN=cluster->GetNumOfDigits();
1932               cnew.fMultiplicity=0;
1933               cnew.fMultiplicityN=nstripsN;
1934               /*
1935               if (nstripsN > 100) {
1936                 printf("multiplicity > 100 - increase the dimension of the arrays %d\n",nstripsN);
1937                 nstripsN=100;
1938               }
1939               Int_t k;
1940               for(k=0;k<nstripsN;k++) {
1941                 // check if 'clusterP->GetDigitStripNo(i)' returns 
1942                 // the digit index and not smth else
1943                 cnew.fIndexMapN[k] = cluster->GetDigitStripNo(k); 
1944               } 
1945               */
1946               cnew.fQErr=sigerr;
1947               //cnew.fProbability=0.75; 
1948               fITS->AddCluster(2,&cnew);
1949               // add the rec point info
1950               AliITSRecPoint rnew;
1951               rnew.SetX(kconv*(x1+x2)/2);
1952               rnew.SetZ(kconv*(z1+z2)/2);
1953               rnew.SetQ(sig);
1954               rnew.SetdEdX(sig/kdEdXtoQ);
1955               rnew.SetSigmaX2( kRMSx* kRMSx);
1956               rnew.SetSigmaZ2( kRMSz* kRMSz);
1957               rnew.SetProbability(0.75);
1958               fITS->AddRecPoint(rnew);
1959               /*
1960               fPointsM->AddLast( (TObject*) 
1961                    new AliITSRecPointSSD((x1+x2)/2 ,0,0,0,(z1+z2)/2,0,sig,sigerr,fLayer,fLad,fDet,cluster, 0.75) );
1962               */
1963
1964               if(debug) cout<<"\n"<<": SINGLE SIDE ImpactPoint Created: X = "
1965                   <<(x1+x2)/2<<" Z = "<<(z1+z2)/2<<";   Pos = "<<pos;
1966  
1967             }
1968         }
1969     }
1970
1971
1972 }
1973
1974 //_______________________________________________________________________
1975
1976 Bool_t AliITSClusterFinderSSD::GetCrossing (Float_t &P, Float_t &N) 
1977
1978
1979    Float_t Dx = fSegmentation->Dx();
1980    Float_t Dz = fSegmentation->Dz();
1981
1982    //P==X
1983    //N==Z
1984
1985     Float_t dx = 0.1;
1986     P *= fPitch;
1987     N *= fPitch; 
1988     
1989     //P = (kZ * fTan + N + P)/2.0;         // x coordinate
1990     //N = kZ - (P-N)/fTan;                 // z coordinate
1991     
1992     if(fTanP + fTanN == 0) return kFALSE;
1993
1994     N = (N - P + fTanN * Dz) / (fTanP + fTanN);  // X coordinate
1995     P = P + fTanP * N;                           // Y coordinate
1996
1997     P -= Dx/2;
1998     N -= Dz/2;
1999
2000     //   N = -(N - P + kZ/2*(fTanP + fTanN))/(fTanP + fTanN);
2001     //  P = (fTanP*(N-kX/2-kZ*fTanN/2) + fTanN*(P-kX/2-kZ*fTanP/2) )/(fTanP + fTanN);
2002     
2003     if ( ( N < -(Dz/2+dx) ) || ( N > (Dz/2+dx) ) ) return kFALSE;
2004     if ( ( P < -(Dx/2+dx) ) || ( P > (Dx/2+dx) ) ) return kFALSE;
2005
2006     return kTRUE;   
2007 }
2008
2009
2010 //_________________________________________________________________________
2011
2012 void AliITSClusterFinderSSD::GetCrossingError(Float_t& dP, Float_t& dN)
2013 {
2014   Float_t dz, dx;
2015   
2016   dz = TMath::Abs(( dP + dN )*fPitch/(fTanP + fTanN) );
2017   dx = fPitch*(TMath::Abs(dP*(1 - fTanP/(fTanP + fTanN))) +
2018                TMath::Abs(dN *fTanP/(fTanP + fTanN) ));
2019   
2020   dN = dz;
2021   dP = dx;
2022 }
2023
2024