Line data Source code
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 118 : ClassImp(AliITSSortTrkl)
34 :
35 : //______________________________________________________________________
36 0 : AliITSSortTrkl::AliITSSortTrkl():TObject(),
37 0 : fkSize(0),
38 0 : fIndex(0),
39 0 : fPairs(NULL),
40 0 : fClustersTmp(NULL),
41 0 : fClusters(NULL),
42 0 : fNoClus(0),
43 0 : fSize(NULL),
44 0 : fMark(),
45 0 : fCut(0.),
46 0 : fCoarseMaxRCut(0.)
47 0 : {
48 : // Default constructor
49 0 : }
50 :
51 : //______________________________________________________________________
52 0 : AliITSSortTrkl::AliITSSortTrkl(Int_t n, Double_t cut):TObject(),
53 0 : fkSize(n),
54 0 : fIndex(0),
55 0 : fPairs(),
56 0 : fClustersTmp(NULL),
57 0 : fClusters(NULL),
58 0 : fNoClus(0),
59 0 : fSize(NULL),
60 0 : fMark(n),
61 0 : fCut(cut),
62 0 : fCoarseMaxRCut(0.) {
63 : // Standard constructor
64 0 : fPairs = new AliITSTracklPairs* [fkSize];
65 0 : fClustersTmp = new Int_t* [fkSize-1];
66 0 : for(Int_t i=0; i<fkSize-1; i++) fClustersTmp[i]=NULL;
67 0 : }
68 :
69 : //______________________________________________________________________
70 0 : AliITSSortTrkl::AliITSSortTrkl(TClonesArray &tclo, Int_t n, Double_t cut, Double_t rcut):TObject(),
71 0 : fkSize(n),
72 0 : fIndex(0),
73 0 : fPairs(),
74 0 : fClustersTmp(NULL),
75 0 : fClusters(NULL),
76 0 : fNoClus(0),
77 0 : fSize(NULL),
78 0 : fMark(n),
79 0 : fCut(cut),
80 0 : fCoarseMaxRCut(rcut){
81 : // Constructor based on a TClonesArray of AliStrLine
82 0 : Int_t numtrack=tclo.GetEntriesFast();
83 0 : if(numtrack<2){
84 0 : AliFatal(Form("Insufficient number of tracks (%d )",numtrack));
85 0 : return;
86 : }
87 0 : TObject *member = tclo[0];
88 0 : TString str(member->ClassName());
89 0 : if(!(str.Contains("AliStrLine"))){
90 0 : AliFatal(Form("Wrong type of class in input TClonesArray (%s )",str.Data()));
91 0 : return;
92 : }
93 0 : Int_t siz = numtrack*(numtrack-1)/2;
94 0 : if(fkSize < siz){
95 0 : AliError(Form("fkSize is too small. It is %d and it should be %d",fkSize,siz));
96 : }
97 0 : fPairs = new AliITSTracklPairs* [fkSize];
98 0 : fClustersTmp = new Int_t* [fkSize-1];
99 0 : for(Int_t i=0; i<fkSize-1; i++) fClustersTmp[i]=NULL;
100 0 : for(Int_t i=0;i<numtrack-1;i++){
101 0 : AliStrLine *one = (AliStrLine*)tclo[i];
102 0 : for(Int_t j=i+1;j<numtrack;j++){
103 0 : AliStrLine *two = (AliStrLine*)tclo[j];
104 0 : Double_t point[3];
105 0 : one->Cross(two,point);
106 0 : Double_t dca = one->GetDCA(two);
107 0 : if(dca>fCut)continue;
108 0 : Double_t rad=TMath::Sqrt(point[0]*point[0]+point[1]*point[1]);
109 0 : if(rad>fCoarseMaxRCut)continue;
110 0 : AddPairs(i,j,dca,point);
111 0 : }
112 : }
113 0 : }
114 :
115 : //______________________________________________________________________
116 0 : AliITSSortTrkl::~AliITSSortTrkl(){
117 : // Destructor
118 0 : if(fPairs){
119 0 : for(Int_t i=0;i<fIndex;i++)delete fPairs[i];
120 0 : delete [] fPairs;
121 : }
122 0 : DeleteClustersTmp();
123 0 : if(fClusters){
124 0 : for(Int_t i=0;i<fNoClus;i++)delete [] fClusters[i];
125 0 : delete [] fClusters;
126 0 : delete [] fSize;
127 : }
128 0 : }
129 :
130 : //______________________________________________________________________
131 : void AliITSSortTrkl::DeleteClustersTmp(){
132 : // fClustersTmp is deleted
133 0 : if(fClustersTmp){
134 0 : for(Int_t i=0;i<fIndex-1;i++){
135 0 : if(fClustersTmp[i]){
136 0 : delete []fClustersTmp[i];
137 : }
138 : }
139 0 : delete [] fClustersTmp;
140 0 : fClustersTmp = NULL;
141 0 : }
142 0 : }
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 0 : 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 0 : fIndex++;
150 0 : return fIndex-1;
151 0 : }
152 :
153 : //______________________________________________________________________
154 : Int_t AliITSSortTrkl::FindClusters(){
155 : // find clusters
156 0 : if(fIndex<2){
157 0 : AliWarning(Form("fIndex = %d",fIndex));
158 0 : fNoClus = 0;
159 0 : return fNoClus;
160 : }
161 0 : fMark.ResetAllBits();
162 0 : PrepareClustersTmp();
163 : // cout<<"AliITSSortTrkl::FindClusters fkSize "<<fkSize<<"; fIndex "<<fIndex<<endl;
164 0 : for(Int_t i=0; i<fIndex-1;i++){
165 0 : Int_t *v=fClustersTmp[i];
166 0 : AliDebug(25,Form("Starting clustering for pair number %d",i));
167 0 : Clustering(i,v);
168 0 : if(v[0]>0){
169 0 : AliDebug(25,Form("Clusters starting from pair %d : %d",i,v[0]));
170 0 : Int_t dim;
171 0 : 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 0 : 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 0 : for(Int_t j=0;j<fIndex;j++){
196 0 : if(fMark.TestBitNumber(j))continue;
197 0 : for(Int_t k=0;k<dim;k++){
198 0 : if(fPairs[j]->HasTrack(arr[k])){
199 0 : fMark.SetBitNumber(j);
200 : // cout<<"Marked pair "<<j<<" because has tracklet "<<arr[k]<<endl;
201 0 : k=dim;
202 0 : }
203 : }
204 0 : }
205 :
206 0 : delete []arr;
207 0 : }
208 : }
209 : // The following method builds the array fClusters
210 0 : Cleanup();
211 0 : return fNoClus;
212 0 : }
213 :
214 : //______________________________________________________________________
215 : void AliITSSortTrkl::Cleanup(){
216 : // empty arrays are eliminated, the others are sorted according
217 : // to cluster multiplicity
218 0 : AliDebug(25,Form("fIndex = %d",fIndex));
219 0 : Int_t *siz = new Int_t[fIndex-1];
220 0 : Int_t *index = new Int_t[fIndex-1];
221 0 : fNoClus=0;
222 0 : for(Int_t i=0; i<fIndex-1;i++){
223 0 : Int_t *v=fClustersTmp[i];
224 0 : if(v[0]>0)fNoClus++;
225 0 : siz[i]=v[0];
226 : }
227 0 : if(fNoClus == 0){
228 0 : delete []siz;
229 0 : delete [] index;
230 0 : return;
231 : }
232 0 : AliDebug(25,Form("fNoClus = %d",fNoClus));
233 0 : TMath::Sort(fIndex-1,siz,index);
234 0 : fClusters = new Int_t* [fNoClus];
235 0 : fSize = new Int_t [fNoClus];
236 0 : for(Int_t i=0;i<fNoClus;i++){
237 : // cout<<"AliITSSortTrkl::Cleanup: Cluster number "<<i<<"; Index= "<<index[i]<<endl;
238 0 : Int_t curind = index[i];
239 0 : Int_t *v=fClustersTmp[curind];
240 0 : fClusters[i] = new Int_t[v[0]+1];
241 0 : fSize[i]=v[0]+1;
242 : // cout<<"AliITSSortTrkl::Cleanup: Size = "<<fSize[i]<<endl;
243 0 : Int_t *vext=fClusters[i];
244 0 : vext[0]=curind;
245 0 : 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 0 : delete []siz;
255 0 : delete [] index;
256 0 : }
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 0 : AliDebug(25,Form("called with parameters %d %d",index,dim));
263 0 : if(fNoClus <=0){
264 0 : dim = 0;
265 0 : return NULL;
266 : }
267 0 : Int_t dimmax = 2*GetSizeOfCluster(index);
268 0 : AliDebug(25,Form("dimmax = %d",dimmax));
269 0 : if(dimmax<=0){
270 0 : dim = 0;
271 0 : return NULL;
272 : }
273 0 : Int_t *v = GetClusters(index);
274 0 : return FindLabels(v,dimmax,dim);
275 0 : }
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 0 : Int_t *arr = new Int_t [dimmax];
285 : Int_t j=0;
286 0 : for(Int_t i=0; i<dimmax/2; i++){
287 0 : AliITSTracklPairs *pai=GetPairsAt(v[i]);
288 0 : arr[j++]=pai->GetTrack1();
289 : // cout<<"AliITSSortTrkl::FindLabels - i="<<i<<" arr["<<j-1<<"]= "<<arr[j-1]<<endl;
290 0 : arr[j++]=pai->GetTrack2();
291 : // cout<<"AliITSSortTrkl::FindLabels arr["<<j-1<<"]= "<<arr[j-1]<<endl;
292 : }
293 0 : SortAndClean(dimmax,arr,dim);
294 0 : 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 0 : if(numb<=0)return;
305 0 : Int_t* index = new Int_t[numb];
306 0 : TMath::Sort(numb,arr,index,kFALSE);
307 0 : Int_t* tmp = new Int_t[numb];
308 0 : numb2 = 0;
309 0 : tmp[0] = arr[index[0]];
310 0 : for(Int_t i=1;i<numb;i++){
311 0 : if(arr[index[i]] != tmp[numb2]){
312 0 : ++numb2;
313 0 : tmp[numb2]=arr[index[i]];
314 0 : }
315 : }
316 0 : ++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 0 : delete [] index;
330 0 : delete [] tmp;
331 0 : }
332 :
333 : //______________________________________________________________________
334 : void AliITSSortTrkl::PrepareClustersTmp(){
335 : // prepare arrays of clusters
336 0 : for(Int_t i=0; i<fIndex-1;i++){
337 0 : fClustersTmp[i] = new Int_t [fIndex+1];
338 0 : Int_t *v = fClustersTmp[i];
339 0 : v[0]=0;
340 0 : for(Int_t j=1;j<fIndex+1;j++)v[j]=-1;
341 : }
342 0 : }
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 0 : AliDebug(25,Form("Clustering called for pair %d",i));
349 0 : if(fMark.TestBitNumber(i)){
350 0 : AliDebug(25,Form("Leaving Clustering for pair %d - nothing done",i));
351 : return;
352 : }
353 0 : AliITSTracklPairs *p1 = fPairs[i];
354 0 : for(Int_t j=0;j<fIndex;j++){
355 0 : if(j == i) continue;
356 0 : if(fMark.TestBitNumber(j))continue;
357 0 : AliITSTracklPairs *p2 = fPairs[j];
358 0 : Double_t dist = p1->GetDistance(*p2);
359 : // AliInfo(Form(" ******* i %d , j %d . Distance %g ",i,j,dist));
360 0 : if(dist<=fCut){
361 0 : Int_t dimclu=v[0];
362 : Bool_t already = kFALSE;
363 0 : for(Int_t k=1;k<=dimclu;k++){
364 0 : if(v[k]==j)already=kTRUE;
365 : }
366 0 : if(!already){
367 0 : ++dimclu;
368 0 : v[0]=dimclu;
369 0 : fMark.SetBitNumber(i);
370 0 : AliDebug(25,Form("Marked pair %d - and call Clustering for pair %d",i,j));
371 0 : v[dimclu]=j;
372 0 : Clustering(j,v);
373 0 : fMark.SetBitNumber(j);
374 0 : AliDebug(25,Form("Marked pair %d",j));
375 : }
376 0 : }
377 0 : }
378 0 : AliDebug(25,Form("Leaving Clustering for pair %d ",i));
379 0 : }
380 :
381 :
382 :
383 :
384 :
|
