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[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 //______________________________________________________________________
 217 Int_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 //______________________________________________________________________
 234 Int_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 //______________________________________________________________________
 249 void 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 //______________________________________________________________________
 278 void 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 //______________________________________________________________________
 289 void 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