Fixing a little overlap. Adding flags to perform material budget studies (Mario)
[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(),
1ca5c3f5 56fClustersTmp(NULL),
108cb8b9 57fClusters(NULL),
58fNoClus(0),
59fSize(NULL),
60fMark(n),
61fCut(cut),
62fCoarseMaxRCut(0.) {
63 // Standard constructor
64 fPairs = new AliITSTracklPairs* [fkSize];
1ca5c3f5 65 fClustersTmp = new Int_t* [fkSize-1];
66 for(Int_t i=0; i<fkSize-1; i++) fClustersTmp[i]=NULL;
108cb8b9 67}
68
69//______________________________________________________________________
70AliITSSortTrkl::AliITSSortTrkl(TClonesArray &tclo, Int_t n, Double_t cut, Double_t rcut):TObject(),
71fkSize(n),
72fIndex(0),
73fPairs(),
1ca5c3f5 74fClustersTmp(NULL),
108cb8b9 75fClusters(NULL),
76fNoClus(0),
77fSize(NULL),
78fMark(n),
79fCut(cut),
80fCoarseMaxRCut(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];
1ca5c3f5 99 for(Int_t i=0; i<fkSize-1; i++) fClustersTmp[i]=NULL;
108cb8b9 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//______________________________________________________________________
116AliITSSortTrkl::~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){
6fefa5ce 124 for(Int_t i=0;i<fNoClus;i++)delete [] fClusters[i];
125 delete [] fClusters;
126 delete [] fSize;
108cb8b9 127 }
128}
129
130//______________________________________________________________________
131void AliITSSortTrkl::DeleteClustersTmp(){
132 // fClustersTmp is deleted
133 if(fClustersTmp){
1ca5c3f5 134 for(Int_t i=0;i<fIndex-1;i++){
135 if(fClustersTmp[i]){
136 delete []fClustersTmp[i];
137 }
138 }
6fefa5ce 139 delete [] fClustersTmp;
108cb8b9 140 fClustersTmp = NULL;
141 }
142}
143
144//______________________________________________________________________
145Int_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//______________________________________________________________________
154Int_t AliITSSortTrkl::FindClusters(){
155 // find clusters
1ca5c3f5 156 if(fIndex<2){
157 AliWarning(Form("fIndex = %d",fIndex));
158 fNoClus = 0;
159 return fNoClus;
160 }
108cb8b9 161 fMark.ResetAllBits();
162 PrepareClustersTmp();
1ca5c3f5 163 // cout<<"AliITSSortTrkl::FindClusters fkSize "<<fkSize<<"; fIndex "<<fIndex<<endl;
108cb8b9 164 for(Int_t i=0; i<fIndex-1;i++){
165 Int_t *v=fClustersTmp[i];
1ca5c3f5 166 AliDebug(25,Form("Starting clustering for pair number %d",i));
167 Clustering(i,v);
108cb8b9 168 if(v[0]>0){
1ca5c3f5 169 AliDebug(25,Form("Clusters starting from pair %d : %d",i,v[0]));
108cb8b9 170 Int_t dim;
171 v[v[0]+1]=i;
1ca5c3f5 172 // arr will contain the labels of the tracks associated to
173 // the cluster of pairs starting from pair i.
108cb8b9 174 Int_t *arr=FindLabels(v+1,2*(v[0]+1),dim);
1ca5c3f5 175
176 /*
177 cout<<"AliITSSortTrkl::FindClusters: Pairs involved \n";
178 for(Int_t j=1;j<=v[0]+1;j++){
108cb8b9 179 cout<<v[j]<<" ";
180 if(j%10==0)cout<<endl;
181 }
182 cout<<endl;
1ca5c3f5 183 cout<<"AliITSSortTrkl::FindClusters: Tracklets involved\n";
108cb8b9 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;
1ca5c3f5 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
108cb8b9 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);
1ca5c3f5 200 // cout<<"Marked pair "<<j<<" because has tracklet "<<arr[k]<<endl;
108cb8b9 201 k=dim;
202 }
203 }
204 }
1ca5c3f5 205
108cb8b9 206 delete []arr;
207 }
208 }
1ca5c3f5 209 // The following method builds the array fClusters
108cb8b9 210 Cleanup();
211 return fNoClus;
212}
213
214//______________________________________________________________________
215void AliITSSortTrkl::Cleanup(){
216 // empty arrays are eliminated, the others are sorted according
217 // to cluster multiplicity
1ca5c3f5 218 AliDebug(25,Form("fIndex = %d",fIndex));
9e817851 219 Int_t *siz = new Int_t[fIndex-1];
220 Int_t *index = new Int_t[fIndex-1];
108cb8b9 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 }
9e817851 227 if(fNoClus == 0){
228 delete []siz;
229 delete [] index;
230 return;
231 }
1ca5c3f5 232 AliDebug(25,Form("fNoClus = %d",fNoClus));
108cb8b9 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++){
1ca5c3f5 237 // cout<<"AliITSSortTrkl::Cleanup: Cluster number "<<i<<"; Index= "<<index[i]<<endl;
108cb8b9 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;
1ca5c3f5 242 // cout<<"AliITSSortTrkl::Cleanup: Size = "<<fSize[i]<<endl;
108cb8b9 243 Int_t *vext=fClusters[i];
244 vext[0]=curind;
245 for(Int_t j=1;j<fSize[i];j++)vext[j]=v[j];
1ca5c3f5 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 */
108cb8b9 253 }
9e817851 254 delete []siz;
255 delete [] index;
108cb8b9 256}
257
258//______________________________________________________________________
259Int_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
1ca5c3f5 262 AliDebug(25,Form("called with parameters %d %d",index,dim));
108cb8b9 263 if(fNoClus <=0){
264 dim = 0;
265 return NULL;
266 }
267 Int_t dimmax = 2*GetSizeOfCluster(index);
1ca5c3f5 268 AliDebug(25,Form("dimmax = %d",dimmax));
108cb8b9 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//______________________________________________________________________
278Int_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
1ca5c3f5 282
283 // cout<<"AliITSSortTrkl::Findlabels parameters "<<v<<" dimmax = "<<dimmax<<" dim "<<dim<<endl;
108cb8b9 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();
1ca5c3f5 289 // cout<<"AliITSSortTrkl::FindLabels - i="<<i<<" arr["<<j-1<<"]= "<<arr[j-1]<<endl;
108cb8b9 290 arr[j++]=pai->GetTrack2();
1ca5c3f5 291 // cout<<"AliITSSortTrkl::FindLabels arr["<<j-1<<"]= "<<arr[j-1]<<endl;
108cb8b9 292 }
293 SortAndClean(dimmax,arr,dim);
294 return arr;
295}
296//______________________________________________________________________
297void 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.
1ca5c3f5 302
303 // cout<<"AliITSSortTrkl::SortAndClean - Parameters: numb= "<<numb<<" numb2= "<<numb2<<endl;
108cb8b9 304 if(numb<=0)return;
9e817851 305 Int_t* index = new Int_t[numb];
108cb8b9 306 TMath::Sort(numb,arr,index,kFALSE);
9e817851 307 Int_t* tmp = new Int_t[numb];
108cb8b9 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;
1ca5c3f5 317 /*
318 cout<<"AliITSSortTrkl::SortAndClean - numb2 = "<<numb2<<endl;
108cb8b9 319 for(Int_t i=0;i<numb;i++){
320 if(i<numb2){
321 arr[i]=tmp[i];
1ca5c3f5 322 cout<<"arr["<<i<<"]= "<<arr[i]<<endl;
108cb8b9 323 }
324 else {
325 arr[i]=0;
326 }
327 }
1ca5c3f5 328 */
9e817851 329 delete [] index;
330 delete [] tmp;
108cb8b9 331}
332
333//______________________________________________________________________
334void AliITSSortTrkl::PrepareClustersTmp(){
335 // prepare arrays of clusters
336 for(Int_t i=0; i<fIndex-1;i++){
1ca5c3f5 337 fClustersTmp[i] = new Int_t [fIndex+1];
108cb8b9 338 Int_t *v = fClustersTmp[i];
339 v[0]=0;
1ca5c3f5 340 for(Int_t j=1;j<fIndex+1;j++)v[j]=-1;
108cb8b9 341 }
342}
343
1ca5c3f5 344
108cb8b9 345//______________________________________________________________________
346void AliITSSortTrkl::Clustering(Int_t i,Int_t *v){
347 // recursive method to build up clusters starting from point i
1ca5c3f5 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 }
108cb8b9 353 AliITSTracklPairs *p1 = fPairs[i];
1ca5c3f5 354 for(Int_t j=0;j<fIndex;j++){
355 if(j == i) continue;
108cb8b9 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;
108cb8b9 369 fMark.SetBitNumber(i);
1ca5c3f5 370 AliDebug(25,Form("Marked pair %d - and call Clustering for pair %d",i,j));
108cb8b9 371 v[dimclu]=j;
372 Clustering(j,v);
1ca5c3f5 373 fMark.SetBitNumber(j);
374 AliDebug(25,Form("Marked pair %d",j));
108cb8b9 375 }
376 }
377 }
1ca5c3f5 378 AliDebug(25,Form("Leaving Clustering for pair %d ",i));
108cb8b9 379}
380
381
382
1ca5c3f5 383
384