70a33a47b5ce0a3541f7d21e63e80d410fe0aa1f
[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 //  * The package was revised and changed by Boris Batiounia in the time     //
18 //  * period of March - June 2001                                            //
19 // **************************************************************************//
20 ///////////////////////////////////////////////////////////////////////////////
21 #include <Riostream.h>
22 #include <TArrayI.h>
23 #include "AliRun.h"
24 #include "AliITS.h"
25 #include "AliITSdigitSSD.h"
26 #include "AliITSRawClusterSSD.h"
27 #include "AliITSRecPoint.h"
28 #include "AliITSMapA1.h"
29 #include "AliITSClusterFinderSSD.h"
30 #include "AliITSclusterSSD.h"
31 #include "AliITSpackageSSD.h"
32 #include "AliITSresponseSSD.h"
33 #include "AliITSsegmentationSSD.h"
34 #include "AliITSgeom.h"
35
36 const Bool_t AliITSClusterFinderSSD::fgkSIDEP=kTRUE;
37 const Bool_t AliITSClusterFinderSSD::fgkSIDEN=kFALSE;
38
39 ClassImp(AliITSClusterFinderSSD)
40
41 //____________________________________________________________________
42 //
43 //  Constructor
44 //______________________________________________________________________
45 AliITSClusterFinderSSD::AliITSClusterFinderSSD():
46 AliITSClusterFinder(),
47 fClusterP(0),
48 fNClusterP(0),
49 fClusterN(0),
50 fNClusterN(0),
51 fPackages(0),
52 fNPackages(0),
53 fDigitsIndexP(0),
54 fNDigitsP(0),
55 fDigitsIndexN(0),
56 fNDigitsN(0),
57 fPitch(0.0),
58 fTanP(0.0),
59 fTanN(0.0),
60 fPNsignalRatio(0.0),
61 fSFF(0),
62 fSFB(0){
63     //Default constructor
64 }
65 //______________________________________________________________________
66 AliITSClusterFinderSSD::AliITSClusterFinderSSD(AliITSsegmentation *seg,
67                                                TClonesArray *digits):
68 AliITSClusterFinder(seg,0),
69 fClusterP(0),
70 fNClusterP(0),
71 fClusterN(0),
72 fNClusterN(0),
73 fPackages(0),
74 fNPackages(0),
75 fDigitsIndexP(0),
76 fNDigitsP(0),
77 fDigitsIndexN(0),
78 fNDigitsN(0),
79 fPitch(0.0),
80 fTanP(0.0),
81 fTanN(0.0),
82 fPNsignalRatio(0.0),
83 fSFF(0),
84 fSFB(0){
85     //Standard constructor
86
87     SetDigits(digits);
88     SetMap(new AliITSMapA1(GetSeg(),Digits()));
89     fClusterP     = new TClonesArray ("AliITSclusterSSD",200);    
90     fNClusterP    = 0;
91     fClusterN     = new TClonesArray ("AliITSclusterSSD",200);   
92     fNClusterN    = 0;
93     fPackages     = new TClonesArray ("AliITSpackageSSD",200);    //packages  
94     fNPackages    = 0;
95     fDigitsIndexP = new TArrayI(300);
96     fNDigitsP     = 0;
97     fDigitsIndexN = new TArrayI(300);
98     fNDigitsN     = 0;
99     fPitch        = GetSeg()->Dpx(0);
100     fPNsignalRatio= 7./8.;    // warning: hard-wired number
101 }
102 //______________________________________________________________________}
103 AliITSClusterFinderSSD::AliITSClusterFinderSSD(AliITSsegmentation *seg,
104                                                AliITSresponse *res):
105 AliITSClusterFinder(seg,res),
106 fClusterP(0),
107 fNClusterP(0),
108 fClusterN(0),
109 fNClusterN(0),
110 fPackages(0),
111 fNPackages(0),
112 fDigitsIndexP(0),
113 fNDigitsP(0),
114 fDigitsIndexN(0),
115 fNDigitsN(0),
116 fPitch(0.0),
117 fTanP(0.0),
118 fTanN(0.0),
119 fPNsignalRatio(0.0),
120 fSFF(0),
121 fSFB(0){
122     //Standard constructor
123
124     fClusterP     = new TClonesArray ("AliITSclusterSSD",200);    
125     fNClusterP    = 0;
126     fClusterN     = new TClonesArray ("AliITSclusterSSD",200);   
127     fNClusterN    = 0;
128     fPackages     = new TClonesArray ("AliITSpackageSSD",200);    //packages  
129     fNPackages    = 0;
130     fDigitsIndexP = new TArrayI(300);
131     fNDigitsP     = 0;
132     fDigitsIndexN = new TArrayI(300);
133     fNDigitsN     = 0;
134     fPitch        = GetSeg()->Dpx(0);
135     fPNsignalRatio= 7./8.;    // warning: hard-wired number
136 }
137
138 //______________________________________________________________________
139 AliITSClusterFinderSSD::AliITSClusterFinderSSD(const AliITSClusterFinderSSD &source) : AliITSClusterFinder(source) {
140   // Copy constructor
141   // Copies are not allowed. The method is protected to avoid misuse.
142   Fatal("AliITSClusterFinderSSD","Copy constructor not allowed\n");
143 }
144
145 //______________________________________________________________________
146 AliITSClusterFinderSSD& AliITSClusterFinderSSD::operator=(const AliITSClusterFinderSSD& /* source */){
147   // Assignment operator
148   // Assignment is not allowed. The method is protected to avoid misuse.
149   Fatal("= operator","Assignment operator not allowed\n");
150   return *this;
151 }
152
153 //______________________________________________________________________
154 AliITSClusterFinderSSD::~AliITSClusterFinderSSD(){
155     // Default destructor
156
157     delete fClusterP;
158     delete fClusterN;        
159     delete fPackages;        
160     delete fDigitsIndexP;        
161     delete fDigitsIndexN; 
162 }
163 //______________________________________________________________________
164 void AliITSClusterFinderSSD::InitReconstruction(){
165     // initialization of the cluster finder
166
167     register Int_t i; //iterator
168
169     for (i=0;i<fNClusterP;i++) fClusterP->RemoveAt(i);
170     fNClusterP  =0;
171     for (i=0;i<fNClusterN;i++) fClusterN->RemoveAt(i);
172     fNClusterN=0;
173     for (i=0;i<fNPackages;i++) fPackages->RemoveAt(i);
174     fNPackages = 0;
175     fNDigitsP  = 0;
176     fNDigitsN  = 0;
177     Float_t stereoP,stereoN;
178     GetSeg()->Angles(stereoP,stereoN);
179     CalcStepFactor(stereoP,stereoN);
180     if(GetDebug(1)) cout<<"fSFF = "<<fSFF<<"  fSFB = "<<fSFB<<"\n";
181 }
182 //______________________________________________________________________
183 void AliITSClusterFinderSSD::FindRawClusters(Int_t module){
184     // This function findes out all clusters belonging to one module
185     // 1. Zeroes all space after previous module reconstruction
186     // 2. Finds all neighbouring digits, create clusters
187     // 3. If necesery, resolves for each group of neighbouring digits 
188     //    how many clusters creates it.
189     // 4. Colculate the x and z coordinate  
190     Int_t lay, lad, detect;
191     AliITSgeom *geom = fITS->GetITSgeom();
192
193     SetModule(module);
194     geom->GetModuleId(GetModule(),lay, lad, detect);
195     if ( lay == 6 ) ((AliITSsegmentationSSD*)GetSeg())->SetLayer(6);
196     if ( lay == 5 ) ((AliITSsegmentationSSD*)GetSeg())->SetLayer(5);
197
198     InitReconstruction();  //ad. 1
199     Map()->FillMap();
200     FillDigitsIndex();
201     SortDigits();
202     FindNeighbouringDigits(); //ad. 2
203     //SeparateOverlappedClusters();  //ad. 3
204     ClustersToPackages();  //ad. 4
205     Map()->ClearMap();
206 }
207 //______________________________________________________________________
208 void AliITSClusterFinderSSD::FindNeighbouringDigits(){
209     //If there are any digits on this side, create 1st Cluster,
210     // add to it this digit, and increment number of clusters
211     register Int_t i;
212
213     if ((fNDigitsP>0 )  && (fNDigitsN > 0 )) {
214         Int_t currentstripNo;
215         Int_t *dbuffer = new Int_t [300];   //buffer for strip numbers
216         Int_t dnumber;    //curent number of digits in buffer
217         TArrayI      &lDigitsIndexP = *fDigitsIndexP;
218         TArrayI      &lDigitsIndexN = *fDigitsIndexN;
219         TObjArray    &lDigits       = *(Digits());
220         TClonesArray &lClusterP     = *fClusterP;
221         TClonesArray &lClusterN     = *fClusterN;
222         //process P side 
223         dnumber = 1;
224         dbuffer[0]=lDigitsIndexP[0];
225         //If next digit is a neigh. of previous, adds to last clust. this digit
226         for (i=1; i<fNDigitsP; i++) {
227             //reads new digit
228             currentstripNo = ((AliITSdigitSSD*)lDigits[lDigitsIndexP[i]])->
229                                                             GetStripNumber(); 
230             //check if it is a neighbour of a previous one
231             if((((AliITSdigitSSD*)lDigits[lDigitsIndexP[i-1]])->
232                                                             GetStripNumber()) 
233                ==  (currentstripNo-1) ) dbuffer[dnumber++]=lDigitsIndexP[i];
234             else{
235                 //create a new one side cluster
236                 new(lClusterP[fNClusterP++]) AliITSclusterSSD(dnumber,dbuffer,
237                                                               Digits(),
238                                                               fgkSIDEP); 
239                 dbuffer[0]=lDigitsIndexP[i];
240                 dnumber = 1;
241             } // end if else
242         } // end loop over fNDigitsP
243         new(lClusterP[fNClusterP++]) AliITSclusterSSD(dnumber,dbuffer,
244                                                       Digits(),fgkSIDEP);
245         //process N side 
246         //for comments, see above
247         dnumber = 1;
248         dbuffer[0]=lDigitsIndexN[0];
249         //If next digit is a neigh. of previous, adds to last clust. this digit
250         for (i=1; i<fNDigitsN; i++) { 
251             currentstripNo = ((AliITSdigitSSD*)(lDigits[lDigitsIndexN[i]]))->
252                                                             GetStripNumber();
253             if ( (((AliITSdigitSSD*)lDigits[lDigitsIndexN[i-1]])->
254                                                             GetStripNumber()) 
255                  == (currentstripNo-1) ) dbuffer[dnumber++]=lDigitsIndexN[i];
256             else {
257                 new(lClusterN[fNClusterN++]) AliITSclusterSSD(dnumber,dbuffer,
258                                                         Digits(),
259                                                         fgkSIDEN);
260                 dbuffer[0]=lDigitsIndexN[i];
261                 dnumber = 1;
262             } // end if else
263         } // end loop over fNDigitsN
264         new(lClusterN[fNClusterN++]) AliITSclusterSSD(dnumber,dbuffer,
265                                                    Digits(),fgkSIDEN);
266         delete [] dbuffer;
267
268     } // end condition on  NDigits 
269
270     if (GetDebug(1)) cout<<"\n Found clusters: fNClusterP = "<<fNClusterP
271                    <<"  fNClusterN ="<<fNClusterN<<"\n";
272 }
273 //______________________________________________________________________
274 void AliITSClusterFinderSSD::SeparateOverlappedClusters(){
275     // overlapped clusters separation
276     register Int_t i; //iterator
277     Double_t  factor=0.75;            // How many percent must be lower signal 
278                                      // on the middle one digit
279                                      // from its neighbours
280     Int_t    signal0;              //signal on the strip before the current one
281     Int_t    signal1;              //signal on the current one signal
282     Int_t    signal2;              //signal on the strip after the current one
283     TArrayI *splitlist;              //  List of splits
284     Int_t    numerofsplits=0;        // number of splits
285     Int_t    initPsize = fNClusterP; //initial size of the arrays 
286     Int_t    initNsize = fNClusterN; //we have to keep it because it will grow 
287                                      // in this function and it doasn't make 
288                                      // sense to pass through it again
289     splitlist = new TArrayI(300);
290
291     for (i=0;i<initPsize;i++){
292         if (( ((AliITSclusterSSD*)(*fClusterP)[i])->
293               GetNumOfDigits())==1) continue;
294         if (( ((AliITSclusterSSD*)(*fClusterP)[i])->
295               GetNumOfDigits())==2) continue;
296         Int_t nj=(((AliITSclusterSSD*)(*fClusterP)[i])->GetNumOfDigits()-1);
297         for (Int_t j=1; j<nj; j++){
298             signal1=((AliITSclusterSSD*)(*fClusterP)[i])->GetDigitSignal(j);
299             signal0=((AliITSclusterSSD*)(*fClusterP)[i])->GetDigitSignal(j-1);
300             signal2=((AliITSclusterSSD*)(*fClusterP)[i])->GetDigitSignal(j+1);
301             //if signal is less then factor*signal of its neighbours
302             if (  (signal1<(factor*signal0)) && (signal1<(factor*signal2)) ){
303                 (*splitlist)[numerofsplits++]=j;
304             } // end if
305         } // end loop over number of digits
306         //split this cluster if necessary
307         if(numerofsplits>0) SplitCluster(splitlist,numerofsplits,i,fgkSIDEP);
308         numerofsplits=0;
309         //in signed places (splitlist)
310     } // end loop over clusters on Pside
311
312     for (i=0;i<initNsize;i++) {
313         if (( ((AliITSclusterSSD*)(*fClusterN)[i])->
314               GetNumOfDigits())==1) continue;
315         if (( ((AliITSclusterSSD*)(*fClusterN)[i])->
316               GetNumOfDigits())==2) continue;
317         Int_t nj=(((AliITSclusterSSD*)(*fClusterN)[i])->GetNumOfDigits()-1);
318         for (Int_t j=1; j<nj; j++){
319             signal1=((AliITSclusterSSD*)(*fClusterN)[i])->GetDigitSignal(j);
320             signal0=((AliITSclusterSSD*)(*fClusterN)[i])->GetDigitSignal(j-1);
321             signal2=((AliITSclusterSSD*)(*fClusterN)[i])->GetDigitSignal(j+1);
322             //if signal is less then factor*signal of its neighbours
323             if (  (signal1<(factor*signal0)) && (signal1<(factor*signal2)) ) 
324                 (*splitlist)[numerofsplits++]=j;  
325         } // end loop over number of digits 
326         //split this cluster into more clusters
327         if(numerofsplits>0) SplitCluster(splitlist,numerofsplits,i,fgkSIDEN);
328         numerofsplits=0;
329         //in signed places (splitlist)
330     } // end loop over clusters on Nside
331
332     delete splitlist;
333 }
334 //______________________________________________________________________
335 void AliITSClusterFinderSSD::SplitCluster(TArrayI *list, Int_t nsplits,
336                                           Int_t index, Bool_t side){
337     //This function splits one side cluster into more clusters
338     //number of splits is defined by "nsplits"
339     //Place of splits are defined in the TArray "list"
340     // For further optimisation: Replace this function by two 
341     // specialised ones (each for one side)
342     // save one "if"
343     //For comlete comments see AliITSclusterSSD::SplitCluster
344     register Int_t i; //iterator
345     AliITSclusterSSD* curentcluster;
346     Int_t   *tmpdigits = new Int_t[100];
347     Int_t    nn;
348
349     // side true means P side
350     if (side) {
351         curentcluster =((AliITSclusterSSD*)((*fClusterP)[index])) ;
352         for (i = nsplits; i>0 ;i--) {  
353             nn=curentcluster->SplitCluster((*list)[(i-1)],tmpdigits);
354             new ((*fClusterP)[fNClusterP]) AliITSclusterSSD(nn,tmpdigits,
355                                                             Digits(),side);
356             ( (AliITSclusterSSD*)((*fClusterP)[fNClusterP]) )->
357                                                       SetLeftNeighbour(kTRUE);
358             //if left cluster had neighbour on the right before split 
359             //new should have it too
360             if ( curentcluster->GetRightNeighbour() ) 
361                 ( (AliITSclusterSSD*)((*fClusterP)[fNClusterP]) )->
362                                                      SetRightNeighbour(kTRUE);
363             else curentcluster->SetRightNeighbour(kTRUE); 
364             fNClusterP++;
365         } // end loop over nplits
366     } else {
367         curentcluster =((AliITSclusterSSD*)((*fClusterN)[index]));
368         for (i = nsplits; i>0 ;i--) {  
369             nn=curentcluster->SplitCluster((*list)[(i-1)],tmpdigits);
370             new ((*fClusterN)[fNClusterN]) AliITSclusterSSD(nn,tmpdigits,
371                                                             Digits(),side);
372             ((AliITSclusterSSD*)((*fClusterN)[fNClusterN]))->
373                                                     SetRightNeighbour(kTRUE);
374             if (curentcluster->GetRightNeighbour())
375                 ( (AliITSclusterSSD*)( (*fClusterN)[fNClusterN]) )->
376                                                      SetRightNeighbour(kTRUE);
377             else curentcluster->SetRightNeighbour(kTRUE);      
378             fNClusterN++;
379         } // end loop over nplits
380     } // end if side
381     delete []tmpdigits;
382 }
383 //______________________________________________________________________
384 Int_t AliITSClusterFinderSSD::SortDigitsP(Int_t start, Int_t end){
385     // sort digits on the P side
386     Int_t right;
387     Int_t left;
388
389     if (start != (end - 1) ){
390         left=this->SortDigitsP(start,(start+end)/2);
391         right=this->SortDigitsP((start+end)/2,end);  
392         return (left || right);
393     }else{ 
394         left =  ((AliITSdigitSSD*)((*(Digits()))[(*fDigitsIndexP)[start]]))->
395                                                               GetStripNumber();
396         right= ((AliITSdigitSSD*)((*(Digits()))[(*fDigitsIndexP)[end]]))->
397                                                               GetStripNumber();
398         if( left > right ){
399             Int_t tmp = (*fDigitsIndexP)[start];
400             (*fDigitsIndexP)[start]=(*fDigitsIndexP)[end];
401             (*fDigitsIndexP)[end]=tmp;
402             return 1;
403         }else return 0;
404     } // end if
405 }
406 //______________________________________________________________________
407 Int_t AliITSClusterFinderSSD::SortDigitsN(Int_t start, Int_t end){
408     // sort digits on the N side
409     Int_t right;
410     Int_t left;
411
412     if (start != (end - 1)){
413         left=this->SortDigitsN(start,(start+end)/2);
414         right=this->SortDigitsN((start+end)/2,end);  
415         return (left || right);
416     }else{
417         left =((AliITSdigitSSD*)((*(Digits()))[(*fDigitsIndexN)[start]]))->
418             GetStripNumber();
419         right=((AliITSdigitSSD*)((*(Digits()))[(*fDigitsIndexN)[end]]))->
420             GetStripNumber();
421         if ( left > right ){
422             Int_t tmp = (*fDigitsIndexN)[start];
423             (*fDigitsIndexN)[start]=(*fDigitsIndexN)[end];
424             (*fDigitsIndexN)[end]=tmp;
425             return 1;
426         }else return 0;
427     } // end if
428 }
429 //______________________________________________________________________
430 void AliITSClusterFinderSSD::FillDigitsIndex(){
431     //Fill the indexes of the clusters belonging to a given ITS module
432     Int_t pns=0, nns=0;
433     Int_t tmp,bit,k;
434     Int_t noentries;
435     Int_t i;
436
437     noentries = NDigits();
438
439     Int_t* psidx = new Int_t [noentries*sizeof(Int_t)];
440     Int_t* nsidx = new Int_t [noentries*sizeof(Int_t)]; 
441     if (fDigitsIndexP==NULL) fDigitsIndexP = new TArrayI(noentries);
442     if (fDigitsIndexN==NULL) fDigitsIndexN = new TArrayI(noentries);
443
444     AliITSdigitSSD *dig;
445
446     for ( i = 0 ; i< noentries; i++ ) {
447         dig = (AliITSdigitSSD*)GetDigit(i);
448         if(dig->IsSideP()) { 
449             bit=1;
450             tmp=dig->GetStripNumber();
451             // I find this totally unnecessary - it's just a 
452             // CPU consuming double check
453             for( k=0;k<pns;k++){
454                 if (tmp==psidx[k]){
455                     if (GetDebug(1)) cout<<"Such a digit exists \n";
456                     bit=0;
457                 } // end if
458             } // end for k
459             // end comment 
460             if(bit) {
461                 fDigitsIndexP->AddAt(i,fNDigitsP++);
462                 psidx[pns++]=tmp;
463             } // end if bit
464         } else {
465             bit=1;
466             tmp=dig->GetStripNumber();
467             // same as above
468             for( k=0;k<nns;k++){
469                 if (tmp==nsidx[k]){
470                     if (GetDebug(1)) cout<<"Such a digit exists \n";
471                     bit=0;
472                 } // end if
473             } // for k
474             // end comment
475             if (bit) {
476                 fDigitsIndexN->AddAt(i,fNDigitsN++);
477                 nsidx[nns++] =tmp;
478             } // end if bit
479         } // end if
480     } // end for i
481
482     delete [] psidx;
483     delete [] nsidx;
484
485     if(GetDebug(1)) cout<<"Digits: P = "<<fNDigitsP<<" N = "<<fNDigitsN<<endl;
486 }
487 //______________________________________________________________________
488 void AliITSClusterFinderSSD::SortDigits(){
489     // sort digits
490     Int_t i;
491
492     if(fNDigitsP>1) for (i=0;i<fNDigitsP-1;i++)
493         if (SortDigitsP(0,(fNDigitsP-1-i))==0) break;
494     if(fNDigitsN>1) for (i=0;i<fNDigitsN-1;i++)
495         if(SortDigitsN(0,(fNDigitsN-1-i))==0) break;
496 }
497 //______________________________________________________________________
498 void AliITSClusterFinderSSD::FillClIndexArrays(Int_t* arrayP,Int_t *arrayN)
499     const{
500     // fill cluster index array
501     register Int_t i;
502
503     for (i=0; i<fNClusterP;i++) arrayP[i]=i;
504     for (i=0; i<fNClusterN;i++) arrayN[i]=i;
505 }
506 //______________________________________________________________________
507 void AliITSClusterFinderSSD::SortClusters(Int_t* arrayP, Int_t *arrayN){
508     // sort clusters
509     Int_t i;
510
511     if(fNClusterP>1) for (i=0;i<fNClusterP-1;i++)
512         if (SortClustersP(0,(fNClusterP-1),arrayP)==0)  break;
513     if(fNClusterN>1) for (i=0;i<fNClusterN-1;i++)
514         if (SortClustersN(0,(fNClusterN-1),arrayN)==0)  break;
515 }
516 //______________________________________________________________________
517 Int_t AliITSClusterFinderSSD::SortClustersP(Int_t start, Int_t end,
518                                             Int_t *array){
519     //Sort P side clusters
520     Int_t right;
521     Int_t left;
522
523     if (start != (end - 1) ) {
524         left=this->SortClustersP(start,(start+end)/2,array);
525         right=this->SortClustersP((start+end)/2,end,array);  
526         return (left || right);
527     } else {
528         left =((AliITSclusterSSD*)((*fClusterP)[array[start]]))->
529             GetDigitStripNo(0);
530         right=((AliITSclusterSSD*)((*fClusterP)[array[ end ]]))->
531             GetDigitStripNo(0);
532         if(left>right) {
533             Int_t tmp = array[start];
534             array[start]=array[end];
535             array[end]=tmp;
536             return 1;
537         } else return 0;
538     } // end if
539 }
540 //______________________________________________________________________
541 Int_t AliITSClusterFinderSSD::SortClustersN(Int_t start, Int_t end, 
542                                             Int_t *array){
543     //Sort N side clusters
544     Int_t right;
545     Int_t left;
546
547     if (start != (end - 1) ) {
548         left=this->SortClustersN(start,(start+end)/2,array);
549         right=this->SortClustersN((start+end)/2,end,array);  
550         return (left || right);
551     } else {
552         left =((AliITSclusterSSD*)((*fClusterN)[array[start]]))->
553             GetDigitStripNo(0);
554         right=((AliITSclusterSSD*)((*fClusterN)[array[ end ]]))->
555             GetDigitStripNo(0);
556         if( left > right) {
557             Int_t tmp = array[start];
558             array[start]=array[end];
559             array[end]=tmp;
560             return 1;
561         } else return 0;
562     } // end if
563 }
564 //______________________________________________________________________
565 void AliITSClusterFinderSSD::ClustersToPackages(){
566     // fill packages   
567     
568     Int_t *oneSclP = new Int_t[fNClusterP];//I want to have sorted 1 S clusters
569     Int_t *oneSclN = new Int_t[fNClusterN];//I can not sort it in TClonesArray
570                                            //so, I create table of indexes and 
571                                            //sort it
572                                            //I do not use TArrayI on purpose
573                                            //MB: well, that's not true that one
574                                            //cannot sort objs in TClonesArray
575     AliITSclusterSSD *currentP;
576     AliITSclusterSSD *currentN;
577     Int_t j1, j2;    
578
579     //Fills in One Side Clusters Index Array
580     FillClIndexArrays(oneSclP,oneSclN); 
581     //Sorts filled Arrays
582     //SortClusters(oneSclP,oneSclN);                   
583
584     fNPackages=1;      
585     new ((*fPackages)[0]) AliITSpackageSSD(fClusterP,fClusterN);
586
587     //This part was includede by Boris Batiounia in March 2001.
588     // Take all recpoint pairs (x coordinates) in both P and N sides  
589     // to calculate z coordinates of the recpoints
590
591     for (j1=0;j1<fNClusterP;j1++) {  
592         currentP = GetPSideCluster(oneSclP[j1]);
593         Double_t xP = currentP->GetPosition();
594         Double_t signalP = currentP->GetTotalSignal();
595         for (j2=0;j2<fNClusterN;j2++) {  
596             currentN = GetNSideCluster(oneSclN[j2]);
597             Double_t xN = currentN->GetPosition();
598             Double_t signalN = currentN->GetTotalSignal();
599             CreateNewRecPoint(xP,1,xN,1,signalP,signalN,currentP,currentN,
600                            0.75);
601         } // end for j2
602     } // end for j1
603
604     delete [] oneSclP;
605     delete [] oneSclN;
606 }
607 //______________________________________________________________________
608 Bool_t AliITSClusterFinderSSD::CreateNewRecPoint(Double_t P,Double_t dP,
609                                                  Double_t N, Double_t dN,
610                                                  Double_t SigP,Double_t SigN, 
611                                                  AliITSclusterSSD *clusterP,
612                                                  AliITSclusterSSD *clusterN,
613                                                  Stat_t prob){
614     // create the recpoints
615     const Double_t kADCtoKeV = 2.16; 
616     // 50 ADC units -> 30000 e-h pairs; 1e-h pair -> 3.6e-3 KeV;
617     // 1 ADC unit -> (30000/50)*3.6e-3 = 2.16 KeV 
618     const Double_t kconv = 1.0e-4;
619     const Double_t kRMSx = 20.0*kconv; 
620     const Double_t kRMSz = 800.0*kconv;
621     Int_t n=0;
622     Int_t *tr;
623     Int_t ntracks;
624
625     if (GetCrossing(P,N)) {
626         //GetCrossingError(dP,dN);
627         dP = dN = prob = 0.0; // to remove unused variable warning.
628         AliITSRawClusterSSD cnew;
629         Int_t nstripsP=clusterP->GetNumOfDigits();
630         Int_t nstripsN=clusterN->GetNumOfDigits();
631         Double_t signal = 0;
632         Double_t dedx = 0;
633         if(SigP>SigN) {
634             signal = SigP;
635             dedx = SigP*kADCtoKeV;
636         }else{
637             signal = SigN;
638             dedx = SigN*kADCtoKeV;
639         } // end if SigP>SigN
640         tr = (Int_t*) clusterP->GetTracks(n);
641         ntracks = clusterP->GetNTracks();
642         cnew.SetSignalP(SigP);
643         cnew.SetSignalN(SigN);
644         cnew.SetMultiplicity(nstripsP);
645         cnew.SetMultN(nstripsN);
646         cnew.SetQErr(TMath::Abs(SigP-SigN));
647         cnew.SetNTrack(ntracks);
648         fITS->AddCluster(2,&cnew);
649         AliITSRecPoint rnew;
650         rnew.SetX(P*kconv);
651         rnew.SetZ(N*kconv);
652         rnew.SetQ(signal);
653         rnew.SetdEdX(dedx);
654         rnew.SetSigmaX2( kRMSx* kRMSx); 
655         rnew.SetSigmaZ2( kRMSz* kRMSz);
656         rnew.fTracks[0]=tr[0];
657         rnew.fTracks[1]=tr[1];
658         rnew.fTracks[2]=tr[2];
659         fITS->AddRecPoint(rnew);
660         return kTRUE;
661     } // end if
662     return kFALSE;  
663 }
664 //______________________________________________________________________
665 void  AliITSClusterFinderSSD::CalcStepFactor(Double_t Psteo, Double_t Nsteo){
666     // calculate the step factor for matching clusters
667     // 95 is the pitch, 4000 - dimension along z ?
668     Double_t dz=GetSeg()->Dz();
669
670     fSFF = ( (Int_t)  (Psteo*dz/fPitch ) );// +1;
671     fSFB = ( (Int_t)  (Nsteo*dz/fPitch ) );// +1;
672 }
673 //______________________________________________________________________
674 AliITSclusterSSD* AliITSClusterFinderSSD::GetPSideCluster(Int_t idx){
675     // get P side clusters
676
677     if((idx<0)||(idx>=fNClusterP)){
678         Info("GetPSideCluster","0<index=%d<=%d out of range",idx,fNClusterP);
679         return 0;
680     }else{
681         return (AliITSclusterSSD*)((*fClusterP)[idx]);
682     } // end if
683 }
684 //______________________________________________________________________
685 AliITSclusterSSD* AliITSClusterFinderSSD::GetNSideCluster(Int_t idx){
686     // get N side clusters
687
688     if((idx<0)||(idx>=fNClusterN)){
689         Info("GetNSideCluster","0<index=%d >= %d out of range",idx,fNClusterN);
690         return 0;
691     }else{
692         return (AliITSclusterSSD*)((*fClusterN)[idx]);
693     } // end if
694 }
695 //______________________________________________________________________
696 Bool_t AliITSClusterFinderSSD::GetCrossing (Double_t &P, Double_t &N){ 
697     // get crossing
698     // This function was rivised and changed by Boris Batiounia in March 2001
699     Double_t dx = GetSeg()->Dx(); // detector size in x direction, microns
700     Double_t dz = GetSeg()->Dz(); // detector size in z direction, microns
701     Double_t xL; // x local coordinate
702     Double_t zL; // z local coordinate
703     Double_t x;  // x = xL + dx/2
704     Double_t z;  // z = zL + dz/2
705     Double_t xP; // x coordinate in the P side from the first P strip
706     Double_t xN; // x coordinate in the N side from the first N strip
707     Float_t stereoP,stereoN;
708
709     GetSeg()->Angles(stereoP,stereoN);
710     fTanP=TMath::Tan(stereoP);
711     fTanN=TMath::Tan(stereoN);
712     Double_t kP = fTanP; // Tangent of 0.0075 mrad
713     Double_t kN = fTanN; // Tangent of 0.0275 mrad
714     P *= fPitch;
715     N *= fPitch; 
716
717     xP = N;      // change the mistake for the P/N
718     xN = P;      // coordinates correspondence in this function
719
720     x = xP + kP*(dz*kN-xP+xN)/(kP+kN);
721     z = (dz*kN-xP+xN)/(kP+kN); 
722     xL = x - dx/2;
723     zL = z - dz/2;
724     P = xL;
725     N = zL;  
726
727     if(TMath::Abs(xL) > dx/2 || TMath::Abs(zL) > dz/2) return kFALSE;
728     
729     // Check that xL and zL are inside the detector for the 
730     // correspondent xP and xN coordinates
731
732     return kTRUE;   
733 }
734 //______________________________________________________________________
735 void AliITSClusterFinderSSD::GetCrossingError(Double_t& dP, Double_t& dN){
736     // get crossing error
737     Double_t dz, dx;
738
739     dz = TMath::Abs(( dP + dN )*fPitch/(fTanP + fTanN) );
740     dx = fPitch*(TMath::Abs(dP*(1 - fTanP/(fTanP + fTanN))) +
741                  TMath::Abs(dN *fTanP/(fTanP + fTanN) ));
742     dN = dz;
743     dP = dx;
744 }