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