New classes for finding multiple vertices (in case of pile-up). They will be used...
[u/mrichter/AliRoot.git] / ITS / AliITSSortTrkl.cxx
CommitLineData
108cb8b9 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
33ClassImp(AliITSSortTrkl)
34
35//______________________________________________________________________
36AliITSSortTrkl::AliITSSortTrkl():TObject(),
37fkSize(0),
38fIndex(0),
39fPairs(NULL),
40fClustersTmp(NULL),
41fClusters(NULL),
42fNoClus(0),
43fSize(NULL),
44fMark(),
45fCut(0.),
46fCoarseMaxRCut(0.)
47{
48 // Default constructor
49}
50
51//______________________________________________________________________
52AliITSSortTrkl::AliITSSortTrkl(Int_t n, Double_t cut):TObject(),
53fkSize(n),
54fIndex(0),
55fPairs(),
56fClustersTmp(),
57fClusters(NULL),
58fNoClus(0),
59fSize(NULL),
60fMark(n),
61fCut(cut),
62fCoarseMaxRCut(0.) {
63 // Standard constructor
64 fPairs = new AliITSTracklPairs* [fkSize];
65 fClustersTmp = new Int_t* [fkSize-1];
66}
67
68//______________________________________________________________________
69AliITSSortTrkl::AliITSSortTrkl(TClonesArray &tclo, Int_t n, Double_t cut, Double_t rcut):TObject(),
70fkSize(n),
71fIndex(0),
72fPairs(),
73fClustersTmp(),
74fClusters(NULL),
75fNoClus(0),
76fSize(NULL),
77fMark(n),
78fCut(cut),
79fCoarseMaxRCut(rcut){
80 // Constructor based on a TClonesArray of AliStrLine
81 Int_t numtrack=tclo.GetEntriesFast();
82 if(numtrack<2){
83 AliFatal(Form("Insufficient number of tracks (%d )",numtrack));
84 return;
85 }
86 TObject *member = tclo[0];
87 TString str(member->ClassName());
88 if(!(str.Contains("AliStrLine"))){
89 AliFatal(Form("Wrong type of class in input TClonesArray (%s )",str.Data()));
90 return;
91 }
92 Int_t siz = numtrack*(numtrack-1)/2;
93 if(fkSize < siz){
94 AliError(Form("fkSize is too small. It is %d and it should be %d",fkSize,siz));
95 }
96 fPairs = new AliITSTracklPairs* [fkSize];
97 fClustersTmp = new Int_t* [fkSize-1];
98 for(Int_t i=0;i<numtrack-1;i++){
99 AliStrLine *one = (AliStrLine*)tclo[i];
100 for(Int_t j=i+1;j<numtrack;j++){
101 AliStrLine *two = (AliStrLine*)tclo[j];
102 Double_t point[3];
103 one->Cross(two,point);
104 Double_t dca = one->GetDCA(two);
105 if(dca>fCut)continue;
106 Double_t rad=TMath::Sqrt(point[0]*point[0]+point[1]*point[1]);
107 if(rad>fCoarseMaxRCut)continue;
108 AddPairs(i,j,dca,point);
109 }
110 }
111}
112
113//______________________________________________________________________
114AliITSSortTrkl::~AliITSSortTrkl(){
115 // Destructor
116 if(fPairs){
117 for(Int_t i=0;i<fIndex;i++)delete fPairs[i];
118 delete [] fPairs;
119 }
120 DeleteClustersTmp();
121 if(fClusters){
122 for(Int_t i=0;i<fNoClus;i++)delete []fClusters[i];
123 delete fClusters;
124 delete []fSize;
125 }
126}
127
128//______________________________________________________________________
129void AliITSSortTrkl::DeleteClustersTmp(){
130 // fClustersTmp is deleted
131 if(fClustersTmp){
132 for(Int_t i=0;i<fIndex-1;i++)delete []fClustersTmp[i];
133 delete fClustersTmp;
134 fClustersTmp = NULL;
135 }
136}
137
138//______________________________________________________________________
139Int_t AliITSSortTrkl::AddPairs(Int_t t1, Int_t t2, Double_t dca, Double_t *coo){
140 // Add a tracklet pair at current position
141 fPairs[fIndex] = new AliITSTracklPairs(t1,t2,dca,coo);
142 // cout<<"Pair "<<fIndex<<" Tracklets "<<t1<<" and "<<t2<<". DCA= "<<dca<<" Crossing "<<coo[0]<<" "<<coo[1]<<" "<<coo[2]<<endl;
143 fIndex++;
144 return fIndex-1;
145}
146
147//______________________________________________________________________
148Int_t AliITSSortTrkl::FindClusters(){
149 // find clusters
150 fMark.ResetAllBits();
151 PrepareClustersTmp();
152 for(Int_t i=0; i<fIndex-1;i++){
153 Int_t *v=fClustersTmp[i];
154 Clustering(i,v);
155 if(v[0]>0){
156 AliInfo(Form("Clusters starting from pair %d : %d",i,v[0]));
157 Int_t dim;
158 v[v[0]+1]=i;
159 Int_t *arr=FindLabels(v+1,2*(v[0]+1),dim);
160 cout<<"Pairs involved \n";
161 for(Int_t j=1;j<=v[0];j++){
162 cout<<v[j]<<" ";
163 if(j%10==0)cout<<endl;
164 }
165 cout<<endl;
166 cout<<"Tracklets involved\n";
167 for(Int_t j=0;j<dim;j++){
168 cout<<arr[j]<<" ";
169 if(j%10==0 && j>0)cout<<endl;
170 }
171 cout<<endl;
172 for(Int_t j=0;j<fIndex;j++){
173 if(fMark.TestBitNumber(j))continue;
174 for(Int_t k=0;k<dim;k++){
175 if(fPairs[j]->HasTrack(arr[k])){
176 fMark.SetBitNumber(j);
177 // cout<<"Marked pair "<<j<<" because has tracklet "<<arr[k]<<endl;
178 k=dim;
179 }
180 }
181 }
182 delete []arr;
183 }
184 }
185 Cleanup();
186 return fNoClus;
187}
188
189//______________________________________________________________________
190void AliITSSortTrkl::Cleanup(){
191 // empty arrays are eliminated, the others are sorted according
192 // to cluster multiplicity
193 Int_t siz[fIndex-1];
194 Int_t index[fIndex-1];
195 fNoClus=0;
196 for(Int_t i=0; i<fIndex-1;i++){
197 Int_t *v=fClustersTmp[i];
198 if(v[0]>0)fNoClus++;
199 siz[i]=v[0];
200 }
201 if(fNoClus == 0)return;
202 TMath::Sort(fIndex-1,siz,index);
203 fClusters = new Int_t* [fNoClus];
204 fSize = new Int_t [fNoClus];
205 for(Int_t i=0;i<fNoClus;i++){
206 Int_t curind = index[i];
207 Int_t *v=fClustersTmp[curind];
208 fClusters[i] = new Int_t[v[0]+1];
209 fSize[i]=v[0]+1;
210 Int_t *vext=fClusters[i];
211 vext[0]=curind;
212 for(Int_t j=1;j<fSize[i];j++)vext[j]=v[j];
213 }
214}
215
216//______________________________________________________________________
217Int_t* AliITSSortTrkl::GetTrackletsLab(Int_t index, Int_t& dim) const {
218 // Returns the tracklet labels corresponding to cluster index
219 // Calling code must take care of memory deallocation
220 if(fNoClus <=0){
221 dim = 0;
222 return NULL;
223 }
224 Int_t dimmax = 2*GetSizeOfCluster(index);
225 if(dimmax<=0){
226 dim = 0;
227 return NULL;
228 }
229 Int_t *v = GetClusters(index);
230 return FindLabels(v,dimmax,dim);
231}
232
233//______________________________________________________________________
234Int_t* AliITSSortTrkl::FindLabels(Int_t *v, Int_t dimmax, Int_t& dim) const {
235 // Returns the tracklet labels corresponding to te list of pairs
236 // contained in v.
237 // Calling code must take care of memory deallocation
238 Int_t *arr = new Int_t [dimmax];
239 Int_t j=0;
240 for(Int_t i=0; i<dimmax/2; i++){
241 AliITSTracklPairs *pai=GetPairsAt(v[i]);
242 arr[j++]=pai->GetTrack1();
243 arr[j++]=pai->GetTrack2();
244 }
245 SortAndClean(dimmax,arr,dim);
246 return arr;
247}
248//______________________________________________________________________
249void AliITSSortTrkl::SortAndClean(Int_t numb, Int_t *arr, Int_t& numb2){
250 // Array arr (with numb elements) is sorted in ascending order.
251 // Then possible reoccurrences
252 // of elements are eliminated. numb2 is the number of remaining elements
253 // after cleanup.
254 if(numb<=0)return;
255 Int_t index[numb];
256 TMath::Sort(numb,arr,index,kFALSE);
257 Int_t tmp[numb];
258 numb2 = 0;
259 tmp[0] = arr[index[0]];
260 for(Int_t i=1;i<numb;i++){
261 if(arr[index[i]] != tmp[numb2]){
262 ++numb2;
263 tmp[numb2]=arr[index[i]];
264 }
265 }
266 ++numb2;
267 for(Int_t i=0;i<numb;i++){
268 if(i<numb2){
269 arr[i]=tmp[i];
270 }
271 else {
272 arr[i]=0;
273 }
274 }
275}
276
277//______________________________________________________________________
278void AliITSSortTrkl::PrepareClustersTmp(){
279 // prepare arrays of clusters
280 for(Int_t i=0; i<fIndex-1;i++){
281 fClustersTmp[i] = new Int_t [fIndex-i];
282 Int_t *v = fClustersTmp[i];
283 v[0]=0;
284 for(Int_t j=1;j<fIndex-i;j++)v[j]=-1;
285 }
286}
287
288//______________________________________________________________________
289void AliITSSortTrkl::Clustering(Int_t i,Int_t *v){
290 // recursive method to build up clusters starting from point i
291 if(fMark.TestBitNumber(i))return;
292 AliITSTracklPairs *p1 = fPairs[i];
293 for(Int_t j=i+1;j<fIndex;j++){
294 if(fMark.TestBitNumber(j))continue;
295 AliITSTracklPairs *p2 = fPairs[j];
296 Double_t dist = p1->GetDistance(*p2);
297 // AliInfo(Form(" ******* i %d , j %d . Distance %g ",i,j,dist));
298 if(dist<=fCut){
299 Int_t dimclu=v[0];
300 Bool_t already = kFALSE;
301 for(Int_t k=1;k<=dimclu;k++){
302 if(v[k]==j)already=kTRUE;
303 }
304 if(!already){
305 ++dimclu;
306 v[0]=dimclu;
307 fMark.SetBitNumber(j);
308 fMark.SetBitNumber(i);
309 v[dimclu]=j;
310 Clustering(j,v);
311 }
312 }
313 }
314}
315
316
317