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(),
  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 }
  67
  68 //______________________________________________________________________
  69 AliITSSortTrkl::AliITSSortTrkl(TClonesArray &tclo, Int_t n, Double_t cut, Double_t rcut):TObject(),
  70 fkSize(n),
  71 fIndex(0),
  72 fPairs(),
  73 fClustersTmp(),
  74 fClusters(NULL),
  75 fNoClus(0),
  76 fSize(NULL),
  77 fMark(n),
  78 fCut(cut),
  79 fCoarseMaxRCut(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 //______________________________________________________________________
 114 AliITSSortTrkl::~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 //______________________________________________________________________
 129 void 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 //______________________________________________________________________
 139 Int_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 //______________________________________________________________________
 148 Int_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 //______________________________________________________________________
 190 void AliITSSortTrkl::Cleanup(){
 191   // empty arrays are eliminated, the others are sorted according
 192   // to cluster multiplicity
 193   Int_t *siz = new Int_t[fIndex-1];
 194   Int_t *index = new Int_t[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){
 202     delete []siz;
 203     delete [] index;
 204     return;
 205   }
 206   TMath::Sort(fIndex-1,siz,index);
 207   fClusters = new Int_t* [fNoClus];
 208   fSize = new Int_t [fNoClus];
 209   for(Int_t i=0;i<fNoClus;i++){
 210     Int_t curind = index[i];
 211     Int_t *v=fClustersTmp[curind];
 212     fClusters[i] = new Int_t[v[0]+1];
 213     fSize[i]=v[0]+1;
 214     Int_t *vext=fClusters[i];
 215     vext[0]=curind;
 216     for(Int_t j=1;j<fSize[i];j++)vext[j]=v[j];
 217   }
 218   delete []siz;
 219   delete [] index;
 220 }
 221
 222 //______________________________________________________________________
 223 Int_t* AliITSSortTrkl::GetTrackletsLab(Int_t index, Int_t& dim) const {
 224   // Returns the tracklet labels corresponding to cluster index
 225   // Calling code must take care of memory deallocation
 226   if(fNoClus <=0){
 227     dim = 0;
 228     return NULL;
 229   }
 230   Int_t dimmax = 2*GetSizeOfCluster(index);
 231   if(dimmax<=0){
 232     dim = 0;
 233     return NULL;
 234   }
 235   Int_t *v = GetClusters(index);
 236   return FindLabels(v,dimmax,dim);
 237 }
 238
 239 //______________________________________________________________________
 240 Int_t* AliITSSortTrkl::FindLabels(Int_t *v, Int_t dimmax, Int_t& dim) const {
 241   // Returns the tracklet labels corresponding to te list of pairs
 242   // contained in v.
 243   // Calling code must take care of memory deallocation
 244   Int_t *arr = new Int_t [dimmax];
 245   Int_t j=0;
 246   for(Int_t i=0; i<dimmax/2; i++){
 247     AliITSTracklPairs *pai=GetPairsAt(v[i]);
 248     arr[j++]=pai->GetTrack1();
 249     arr[j++]=pai->GetTrack2();
 250   }
 251   SortAndClean(dimmax,arr,dim);
 252   return arr;
 253 }
 254 //______________________________________________________________________
 255 void AliITSSortTrkl::SortAndClean(Int_t numb, Int_t *arr, Int_t& numb2){
 256   // Array arr (with numb elements) is sorted in ascending order.
 257   // Then possible reoccurrences
 258   // of elements are eliminated. numb2 is the number of remaining elements
 259   // after cleanup.
 260   if(numb<=0)return;
 261   Int_t* index = new Int_t[numb];
 262   TMath::Sort(numb,arr,index,kFALSE);
 263   Int_t* tmp = new Int_t[numb];
 264   numb2 = 0;
 265   tmp[0] = arr[index[0]];
 266   for(Int_t i=1;i<numb;i++){
 267     if(arr[index[i]] != tmp[numb2]){
 268       ++numb2;
 269       tmp[numb2]=arr[index[i]];
 270     }
 271   }
 272   ++numb2;
 273   for(Int_t i=0;i<numb;i++){
 274     if(i<numb2){
 275       arr[i]=tmp[i];
 276     }
 277     else {
 278       arr[i]=0;
 279     }
 280   }
 281   delete [] index;
 282   delete [] tmp;
 283 }
 284
 285 //______________________________________________________________________
 286 void AliITSSortTrkl::PrepareClustersTmp(){
 287   // prepare arrays of clusters
 288   for(Int_t i=0; i<fIndex-1;i++){
 289     fClustersTmp[i] = new Int_t [fIndex-i];
 290     Int_t *v = fClustersTmp[i];
 291     v[0]=0;
 292     for(Int_t j=1;j<fIndex-i;j++)v[j]=-1;
 293   }
 294 }
 295
 296 //______________________________________________________________________
 297 void AliITSSortTrkl::Clustering(Int_t i,Int_t *v){
 298   // recursive method to build up clusters starting from point i
 299   if(fMark.TestBitNumber(i))return;
 300   AliITSTracklPairs *p1 = fPairs[i];
 301   for(Int_t j=i+1;j<fIndex;j++){
 302     if(fMark.TestBitNumber(j))continue;
 303     AliITSTracklPairs *p2 = fPairs[j];
 304     Double_t dist = p1->GetDistance(*p2);
 305     //    AliInfo(Form("  ******* i %d , j %d . Distance %g ",i,j,dist));
 306     if(dist<=fCut){
 307       Int_t dimclu=v[0];
 308       Bool_t already = kFALSE;
 309       for(Int_t k=1;k<=dimclu;k++){
 310         if(v[k]==j)already=kTRUE;
 311       }
 312       if(!already){
 313         ++dimclu;
 314         v[0]=dimclu;
 315         fMark.SetBitNumber(j);
 316         fMark.SetBitNumber(i);
 317         v[dimclu]=j;
 318         Clustering(j,v);
 319       }
 320     }
 321   }
 322 }
 323
 324
 325