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