Correct treatment of tracks with pT<pTmin
[u/mrichter/AliRoot.git] / ITS / AliITSSortTrkl.cxx
1 /**************************************************************************
2  * Copyright(c) 2009-2010, 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 #include<Riostream.h>
16 #include<TClonesArray.h>
17 #include<TMath.h>
18 #include "AliStrLine.h"
19 #include "AliITSSortTrkl.h"
20
21 /* $Id$ */
22
23 ////////////////////////////////////////////////////////////////////////
24 //           Helper class for finding multiple primary vertices       //
25 //           To be used by AliITSVertexer3D                           //
26 //           It is based on the association of pairs of tracklets     //
27 //           obtained by matching reconstructed points onthe first    //
28 //           2 layers of SPD                                          //
29 //           Origin M. Masera masera@to.infn.it                       //
30 ////////////////////////////////////////////////////////////////////////
31
32
33 ClassImp(AliITSSortTrkl)
34
35 //______________________________________________________________________
36 AliITSSortTrkl::AliITSSortTrkl():TObject(),
37 fkSize(0),
38 fIndex(0),
39 fPairs(NULL),
40 fClustersTmp(NULL),
41 fClusters(NULL),
42 fNoClus(0),
43 fSize(NULL),
44 fMark(),
45 fCut(0.),
46 fCoarseMaxRCut(0.)
47 {
48   // Default constructor
49 }
50
51 //______________________________________________________________________
52 AliITSSortTrkl::AliITSSortTrkl(Int_t n, Double_t cut):TObject(),
53 fkSize(n),
54 fIndex(0),
55 fPairs(),
56 fClustersTmp(NULL),
57 fClusters(NULL),
58 fNoClus(0),
59 fSize(NULL),
60 fMark(n),
61 fCut(cut),
62 fCoarseMaxRCut(0.) {
63   // Standard constructor
64   fPairs = new AliITSTracklPairs* [fkSize];
65   fClustersTmp = new Int_t* [fkSize-1];
66   for(Int_t i=0; i<fkSize-1; i++) fClustersTmp[i]=NULL;
67 }
68
69 //______________________________________________________________________
70 AliITSSortTrkl::AliITSSortTrkl(TClonesArray &tclo, Int_t n, Double_t cut, Double_t rcut):TObject(),
71 fkSize(n),
72 fIndex(0),
73 fPairs(),
74 fClustersTmp(NULL),
75 fClusters(NULL),
76 fNoClus(0),
77 fSize(NULL),
78 fMark(n),
79 fCut(cut),
80 fCoarseMaxRCut(rcut){
81   // Constructor based on a TClonesArray of AliStrLine
82   Int_t numtrack=tclo.GetEntriesFast();
83   if(numtrack<2){
84     AliFatal(Form("Insufficient number of tracks (%d )",numtrack));
85     return;
86   }
87   TObject *member = tclo[0];
88   TString str(member->ClassName());
89   if(!(str.Contains("AliStrLine"))){
90     AliFatal(Form("Wrong type of class in input TClonesArray (%s )",str.Data()));
91     return;
92   }
93   Int_t siz = numtrack*(numtrack-1)/2;
94   if(fkSize < siz){
95     AliError(Form("fkSize is too small. It is %d and it should be %d",fkSize,siz)); 
96   }
97   fPairs = new AliITSTracklPairs* [fkSize];
98   fClustersTmp = new Int_t* [fkSize-1];
99   for(Int_t i=0; i<fkSize-1; i++) fClustersTmp[i]=NULL;
100   for(Int_t i=0;i<numtrack-1;i++){
101     AliStrLine *one = (AliStrLine*)tclo[i];
102     for(Int_t j=i+1;j<numtrack;j++){
103       AliStrLine *two = (AliStrLine*)tclo[j];
104       Double_t point[3];
105       one->Cross(two,point);
106       Double_t dca = one->GetDCA(two);
107       if(dca>fCut)continue;
108       Double_t rad=TMath::Sqrt(point[0]*point[0]+point[1]*point[1]);
109       if(rad>fCoarseMaxRCut)continue;
110       AddPairs(i,j,dca,point);
111     }
112   } 
113 }
114
115 //______________________________________________________________________
116 AliITSSortTrkl::~AliITSSortTrkl(){
117   // Destructor
118   if(fPairs){
119     for(Int_t i=0;i<fIndex;i++)delete fPairs[i];
120     delete [] fPairs;
121   }
122   DeleteClustersTmp();
123   if(fClusters){
124     for(Int_t i=0;i<fNoClus;i++)delete [] fClusters[i];
125     delete [] fClusters;
126     delete [] fSize;
127   } 
128 }
129
130 //______________________________________________________________________
131 void AliITSSortTrkl::DeleteClustersTmp(){
132   // fClustersTmp is deleted
133   if(fClustersTmp){
134     for(Int_t i=0;i<fIndex-1;i++){
135       if(fClustersTmp[i]){
136         delete []fClustersTmp[i];
137       }
138     }
139     delete [] fClustersTmp;
140     fClustersTmp = NULL;
141   }
142 }
143
144 //______________________________________________________________________
145 Int_t AliITSSortTrkl::AddPairs(Int_t t1, Int_t t2, Double_t dca, Double_t *coo){
146   // Add a tracklet pair at current position
147   fPairs[fIndex] = new AliITSTracklPairs(t1,t2,dca,coo);
148   //  cout<<"Pair "<<fIndex<<" Tracklets "<<t1<<" and "<<t2<<". DCA= "<<dca<<" Crossing "<<coo[0]<<" "<<coo[1]<<" "<<coo[2]<<endl;
149   fIndex++;
150   return fIndex-1;
151 }
152
153 //______________________________________________________________________
154 Int_t AliITSSortTrkl::FindClusters(){
155   // find clusters
156   if(fIndex<2){
157     AliWarning(Form("fIndex = %d",fIndex));
158     fNoClus = 0;
159     return fNoClus;
160   }
161   fMark.ResetAllBits();
162   PrepareClustersTmp();
163   //  cout<<"AliITSSortTrkl::FindClusters fkSize "<<fkSize<<"; fIndex "<<fIndex<<endl;
164   for(Int_t i=0; i<fIndex-1;i++){
165     Int_t *v=fClustersTmp[i];
166     AliDebug(25,Form("Starting clustering for pair number %d",i));
167     Clustering(i,v);  
168     if(v[0]>0){
169       AliDebug(25,Form("Clusters starting from pair %d : %d",i,v[0]));
170       Int_t dim;
171       v[v[0]+1]=i;
172       // arr will contain the labels of the tracks associated to
173       // the cluster of pairs starting from pair i.
174       Int_t *arr=FindLabels(v+1,2*(v[0]+1),dim);
175
176       /*
177       cout<<"AliITSSortTrkl::FindClusters: Pairs involved \n";
178       for(Int_t j=1;j<=v[0]+1;j++){
179         cout<<v[j]<<" ";
180         if(j%10==0)cout<<endl;
181       }
182       cout<<endl;
183       cout<<"AliITSSortTrkl::FindClusters: Tracklets involved\n";
184       for(Int_t j=0;j<dim;j++){
185         cout<<arr[j]<<" ";
186         if(j%10==0 && j>0)cout<<endl;
187       }
188       cout<<endl;
189       */
190
191       //In the following loop, all the pairs having 
192       // one tracklet already associated to the cluster starting
193       // at pair i are marked, in order to be excluded from further
194       // searches    
195       for(Int_t j=0;j<fIndex;j++){
196         if(fMark.TestBitNumber(j))continue;
197         for(Int_t k=0;k<dim;k++){
198           if(fPairs[j]->HasTrack(arr[k])){
199             fMark.SetBitNumber(j);
200             // cout<<"Marked pair "<<j<<" because has tracklet "<<arr[k]<<endl;
201             k=dim;
202           }
203         }
204       }
205     
206       delete []arr;
207     }
208   }
209   // The following method builds the array fClusters
210   Cleanup();
211   return fNoClus;
212 }
213
214 //______________________________________________________________________
215 void AliITSSortTrkl::Cleanup(){
216   // empty arrays are eliminated, the others are sorted according
217   // to cluster multiplicity
218   AliDebug(25,Form("fIndex = %d",fIndex));
219   Int_t *siz = new Int_t[fIndex-1];
220   Int_t *index = new Int_t[fIndex-1];
221   fNoClus=0;
222   for(Int_t i=0; i<fIndex-1;i++){
223     Int_t *v=fClustersTmp[i];
224     if(v[0]>0)fNoClus++;
225     siz[i]=v[0];
226   }
227   if(fNoClus == 0){
228     delete []siz;
229     delete [] index;
230     return;
231   }
232   AliDebug(25,Form("fNoClus = %d",fNoClus));
233   TMath::Sort(fIndex-1,siz,index);
234   fClusters = new Int_t* [fNoClus];
235   fSize = new Int_t [fNoClus];
236   for(Int_t i=0;i<fNoClus;i++){
237     //    cout<<"AliITSSortTrkl::Cleanup: Cluster number "<<i<<"; Index= "<<index[i]<<endl;
238     Int_t curind = index[i];
239     Int_t *v=fClustersTmp[curind];
240     fClusters[i] = new Int_t[v[0]+1];
241     fSize[i]=v[0]+1;
242     //    cout<<"AliITSSortTrkl::Cleanup: Size = "<<fSize[i]<<endl;
243     Int_t *vext=fClusters[i];
244     vext[0]=curind;
245     for(Int_t j=1;j<fSize[i];j++)vext[j]=v[j];
246     /*
247     for(Int_t j=0;j<fSize[i];j++){
248       cout<<vext[j]<<" ";
249       if(j%10 == 0 && j!=0)cout<<endl;
250     }
251     cout<<endl;
252     */
253   }
254   delete []siz;
255   delete [] index;
256 }
257
258 //______________________________________________________________________
259 Int_t* AliITSSortTrkl::GetTrackletsLab(Int_t index, Int_t& dim) const {
260   // Returns the tracklet labels corresponding to cluster index
261   // Calling code must take care of memory deallocation
262   AliDebug(25,Form("called with parameters %d %d",index,dim));
263   if(fNoClus <=0){
264     dim = 0;
265     return NULL;
266   }
267   Int_t dimmax = 2*GetSizeOfCluster(index);
268   AliDebug(25,Form("dimmax = %d",dimmax));
269   if(dimmax<=0){
270     dim = 0;
271     return NULL;
272   } 
273   Int_t *v = GetClusters(index);
274   return FindLabels(v,dimmax,dim);
275 }
276
277 //______________________________________________________________________
278 Int_t* AliITSSortTrkl::FindLabels(Int_t *v, Int_t dimmax, Int_t& dim) const {
279   // Returns the tracklet labels corresponding to te list of pairs 
280   // contained in v.
281   // Calling code must take care of memory deallocation
282
283   //  cout<<"AliITSSortTrkl::Findlabels parameters "<<v<<" dimmax = "<<dimmax<<" dim "<<dim<<endl;
284   Int_t *arr = new Int_t [dimmax];
285   Int_t j=0;
286   for(Int_t i=0; i<dimmax/2; i++){
287     AliITSTracklPairs *pai=GetPairsAt(v[i]);
288     arr[j++]=pai->GetTrack1();
289     //    cout<<"AliITSSortTrkl::FindLabels - i="<<i<<" arr["<<j-1<<"]= "<<arr[j-1]<<endl;
290     arr[j++]=pai->GetTrack2();
291     //    cout<<"AliITSSortTrkl::FindLabels arr["<<j-1<<"]= "<<arr[j-1]<<endl;
292   }
293   SortAndClean(dimmax,arr,dim);
294   return arr;
295 }
296 //______________________________________________________________________
297 void AliITSSortTrkl::SortAndClean(Int_t numb, Int_t *arr, Int_t& numb2){
298   // Array arr (with numb elements) is sorted in ascending order. 
299   // Then possible reoccurrences
300   // of elements are eliminated. numb2 is the number of remaining elements
301   // after cleanup.
302    
303   //  cout<<"AliITSSortTrkl::SortAndClean - Parameters: numb= "<<numb<<" numb2= "<<numb2<<endl;
304   if(numb<=0)return;
305   Int_t* index = new Int_t[numb];
306   TMath::Sort(numb,arr,index,kFALSE);
307   Int_t* tmp = new Int_t[numb];
308   numb2 = 0;
309   tmp[0] = arr[index[0]];
310   for(Int_t i=1;i<numb;i++){
311     if(arr[index[i]] != tmp[numb2]){
312       ++numb2;
313       tmp[numb2]=arr[index[i]];
314     }
315   }
316   ++numb2;
317   /*
318   cout<<"AliITSSortTrkl::SortAndClean - numb2 = "<<numb2<<endl;
319   for(Int_t i=0;i<numb;i++){
320     if(i<numb2){
321       arr[i]=tmp[i];
322       cout<<"arr["<<i<<"]= "<<arr[i]<<endl;
323     }
324     else {
325       arr[i]=0;
326     }
327   }
328   */
329   delete [] index;
330   delete [] tmp;
331 }
332
333 //______________________________________________________________________
334 void AliITSSortTrkl::PrepareClustersTmp(){
335   // prepare arrays of clusters
336   for(Int_t i=0; i<fIndex-1;i++){
337     fClustersTmp[i] = new Int_t [fIndex+1];
338     Int_t *v = fClustersTmp[i];
339     v[0]=0;
340     for(Int_t j=1;j<fIndex+1;j++)v[j]=-1;
341   }
342 }
343
344
345 //______________________________________________________________________
346 void AliITSSortTrkl::Clustering(Int_t i,Int_t *v){
347   // recursive method to build up clusters starting from point i
348   AliDebug(25,Form("Clustering called for pair %d",i));
349   if(fMark.TestBitNumber(i)){
350     AliDebug(25,Form("Leaving Clustering for pair %d - nothing done",i));
351   return;
352   }
353   AliITSTracklPairs *p1 = fPairs[i];
354   for(Int_t j=0;j<fIndex;j++){
355     if(j == i) continue;
356     if(fMark.TestBitNumber(j))continue;
357     AliITSTracklPairs *p2 = fPairs[j];
358     Double_t dist = p1->GetDistance(*p2);
359     //    AliInfo(Form("  ******* i %d , j %d . Distance %g ",i,j,dist));
360     if(dist<=fCut){
361       Int_t dimclu=v[0];
362       Bool_t already = kFALSE;
363       for(Int_t k=1;k<=dimclu;k++){
364         if(v[k]==j)already=kTRUE;
365       }
366       if(!already){
367         ++dimclu;
368         v[0]=dimclu;
369         fMark.SetBitNumber(i);
370         AliDebug(25,Form("Marked pair %d - and call Clustering for pair %d",i,j));
371         v[dimclu]=j;
372         Clustering(j,v);
373         fMark.SetBitNumber(j);
374         AliDebug(25,Form("Marked pair %d",j));
375       }
376     }
377   }
378   AliDebug(25,Form("Leaving Clustering for pair %d ",i));
379 }
380
381
382
383
384