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