ITS/AliITSClusterFinder.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, 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 ////////////////////////////////////////////////////////////////////////////
  16 //                                                                        //
  17 // Base Class used to find                                                //
  18 // the reconstructed points for ITS                                       //
  19 // See also AliITSClusterFinderSPD, AliITSClusterFinderSDD,               //
  20 // AliITSClusterFinderSDD  AliITSClusterFinderV2                          //
  21 ////////////////////////////////////////////////////////////////////////////
  22
  23 #include "AliRun.h"
  24 #include "AliITSClusterFinder.h"
  25 #include "AliITSRecPoint.h"
  26 #include "AliITSdigit.h"
  27 #include "AliITSDetTypeRec.h"
  28 #include "AliITSMap.h"
  29 #include "AliITSgeomTGeo.h"
  30 #include <TParticle.h>
  31 #include "AliMC.h"
  32
  33 ClassImp(AliITSClusterFinder)
  34
  35 extern AliRun *gAlice;
  36
  37 //----------------------------------------------------------------------
  38 AliITSClusterFinder::AliITSClusterFinder():
  39 TObject(),
  40 fModule(0),
  41 fDigits(0),
  42 fNdigits(0),
  43 fDetTypeRec(0),
  44 fClusters(0),
  45 fMap(0),
  46 fNPeaks(-1),
  47 fNModules(AliITSgeomTGeo::GetNModules()),
  48 fEvent(0){
  49     // default cluster finder
  50     // Input:
  51     //   none.
  52     // Output:
  53     //   none.
  54     // Return:
  55     //   A default constructed AliITSCulsterFinder
  56 }
  57 //----------------------------------------------------------------------
  58 AliITSClusterFinder::AliITSClusterFinder(AliITSDetTypeRec* dettyp):
  59 TObject(),
  60 fModule(0),
  61 fDigits(0),
  62 fNdigits(0),
  63 fDetTypeRec(dettyp),
  64 fClusters(0),
  65 fMap(0),
  66 fNPeaks(-1),
  67 fNModules(AliITSgeomTGeo::GetNModules()),
  68 fEvent(0){
  69     // default cluster finder
  70     // Standard constructor for cluster finder
  71     // Input:
  72     //   AliITSsegmentation *seg  The segmentation class to be used
  73     //   AliITSresponse     *res  The response class to be used
  74     // Output:
  75     //   none.
  76     // Return:
  77     //   A Standard constructed AliITSCulsterFinder
  78
  79 }
  80 //----------------------------------------------------------------------
  81 AliITSClusterFinder::AliITSClusterFinder(AliITSDetTypeRec* dettyp,
  82                                          TClonesArray *digits):
  83 TObject(),
  84 fModule(0),
  85 fDigits(digits),
  86 fNdigits(0),
  87 fDetTypeRec(dettyp),
  88 fClusters(0),
  89 fMap(0),
  90 fNPeaks(-1),
  91 fNModules(AliITSgeomTGeo::GetNModules()),
  92 fEvent(0){
  93     // default cluster finder
  94     // Standard + cluster finder constructor
  95     // Input:
  96     //   AliITSsegmentation *seg  The segmentation class to be used
  97     //   AliITSresponse     *res  The response class to be used
  98     //   TClonesArray    *digits  Array of digits to be used
  99     // Output:
 100     //   none.
 101     // Return:
 102     //   A Standard constructed AliITSCulsterFinder
 103
 104     fNdigits = fDigits->GetEntriesFast();
 105 }
 106
 107 //______________________________________________________________________
 108 AliITSClusterFinder::AliITSClusterFinder(const AliITSClusterFinder &source) : TObject(source),
 109 fModule(source.fModule),
 110 fDigits(),
 111 fNdigits(source.fNdigits),
 112 fDetTypeRec(),
 113 fClusters(),
 114 fMap(),
 115 fNPeaks(source.fNPeaks),
 116 fNModules(source.fNModules),
 117 fEvent(source.fEvent) {
 118   // Copy constructor
 119   // Copies are not allowed. The method is protected to avoid misuse.
 120   AliError("Copy constructor not allowed\n");
 121 }
 122
 123
 124 //______________________________________________________________________
 125 //AliITSClusterFinder& AliITSClusterFinder::operator=(const AliITSClusterFinder& /* source */){
 126   // Assignment operator
 127   // Assignment is not allowed. The method is protected to avoid misuse.
 128 //  Fatal("= operator","Assignment operator not allowed\n");
 129 //  return *this;
 130 //}
 131
 132 //----------------------------------------------------------------------
 133 AliITSClusterFinder::~AliITSClusterFinder(){
 134     // destructor cluster finder
 135     // Input:
 136     //   none.
 137     // Output:
 138     //   none.
 139     // Return:
 140     //   none.
 141
 142     if(fMap) {delete fMap;}
 143     // Zero local pointers. Other classes own these pointers.
 144     fMap          = 0;
 145     fDigits       = 0;
 146     fNdigits      = 0;
 147     fNPeaks       = 0;
 148     fDetTypeRec   = 0;
 149
 150 }
 151 //__________________________________________________________________________
 152 void AliITSClusterFinder::InitGeometry(){
 153  //
 154  // Initialisation of ITS geometry
 155  //
 156   Int_t mmax=AliITSgeomTGeo::GetNModules();
 157   for (Int_t m=0; m<mmax; m++) {
 158     Int_t lay,lad,det; AliITSgeomTGeo::GetModuleId(m,lay,lad,det);
 159     fNdet[m] = (lad-1)*AliITSgeomTGeo::GetNDetectors(lay) + (det-1);
 160     fNlayer[m] = lay-1;
 161   }
 162 }
 163
 164
 165
 166
 167 //______________________________________________________________________
 168 Bool_t AliITSClusterFinder::IsNeighbor(TObjArray *digs,Int_t i,Int_t n[])const{
 169     // Locagical function which checks to see if digit i has a neighbor.
 170     // If so, then it returns kTRUE and its neighbor index j.
 171     // This routine checks if the digits are side by side or one before the
 172     // other. Requires that the array of digits be in proper order.
 173     // Returns kTRUE in the following cases.
 174     //                 ji   0j   if kdiagonal  j0    0i
 175     //                 00   0i   if kdiagonal  0i    j0
 176     // Inputs:
 177     //    TObjArray *digs   Array to search for neighbors in
 178     //    Int_t      i      Index of digit for which we are searching for
 179     //                      a neighbor of.
 180     // Output:
 181     //    Int_t      j[4]   Index of one or more of the digits which is a
 182     //                      neighbor of digit a index i.
 183     // Return:
 184     //    Bool_t            kTRUE if a neighbor was found kFALSE otherwise.
 185     Int_t ix,jx,iz,jz,j;
 186     const Bool_t kdiagonal=kFALSE;
 187     Bool_t nei[4];
 188
 189     // No neighbors found if array empty.
 190     if(digs->GetEntriesFast()<=0) return kFALSE;
 191     // can not be a digit with first element or elements out or range
 192     if(i<=0 || i>=digs->GetEntriesFast()) return kFALSE;
 193
 194     for(j=0;j<4;j++){n[j] = -1;nei[j]=kFALSE;}
 195     ix = ((AliITSdigit*)(digs->At(i)))->GetCoord1();
 196     iz = ((AliITSdigit*)(digs->At(i)))->GetCoord2();
 197     for(j=0;j<i;j++){
 198         jx = ((AliITSdigit*)(digs->At(j)))->GetCoord1();
 199         jz = ((AliITSdigit*)(digs->At(j)))->GetCoord2();
 200         if(jx+1==ix && jz  ==iz){n[0] = j;nei[0] = kTRUE;}
 201         if(jx  ==ix && jz+1==iz){n[1] = j;nei[1] = kTRUE;}
 202         if(jx+1==ix && jz+1==iz){n[2] = j;nei[2] = kTRUE;}
 203         if(jx+1==ix && jz-1==iz){n[3] = j;nei[3] = kTRUE;}
 204     } // end for k
 205     if(nei[0]||nei[1]) return kTRUE;
 206     if(kdiagonal&&(nei[2]||nei[3])) return kTRUE;
 207     // no Neighbors found.
 208     return kFALSE;
 209 }
 210
 211 //______________________________________________________________________
 212 void AliITSClusterFinder::Print(ostream *os) const{
 213     //Standard output format for this class
 214     // Inputs:
 215     //    ostream *os   Output stream
 216     // Output:
 217     //    ostream *os   Output stream
 218     // Return:
 219     //    none.
 220
 221     *os << fModule<<",";
 222     *os << fNdigits<<",";
 223     *os << fNPeaks<<endl;
 224 }
 225 //______________________________________________________________________
 226 void AliITSClusterFinder::Read(istream *is)  {
 227     //Standard input for this class
 228     // Inputs:
 229     //    istream *is   Input stream
 230     // Output:
 231     //    istream *is   Input stream
 232     // Return:
 233     //    none.
 234
 235     *is >> fModule;
 236     *is >> fNdigits;
 237     *is >> fNPeaks;
 238 }
 239 //______________________________________________________________________
 240 ostream &operator<<(ostream &os,AliITSClusterFinder &source){
 241     // Standard output streaming function.
 242     // Inputs:
 243     //    ostream             *os     Output stream
 244     //    AliITSClusterFinder &source Class to be printed
 245     // Output:
 246     //    ostream             *os     Output stream
 247     // Return:
 248     //    none.
 249
 250     source.Print(&os);
 251     return os;
 252 }
 253 //______________________________________________________________________
 254 istream &operator>>(istream &is,AliITSClusterFinder &source){
 255     // Standard output streaming function.
 256     // Inputs:
 257     //    istream              *is      Input stream
 258     //     AliITSClusterFinder &source  Class to be read in.
 259     // Output:
 260     //    istream              *is      Input stream
 261     // Return:
 262     //    none.
 263
 264     source.Read(&is);
 265     return is;
 266 }
 267 //______________________________________________________________________
 268 void AliITSClusterFinder::CheckLabels2(Int_t lab[10]) {
 269   //------------------------------------------------------------
 270   // Tries to find mother's labels
 271   //------------------------------------------------------------
 272   AliRunLoader *rl = AliRunLoader::Instance();
 273   if(!rl) return;
 274   TTree *trK=(TTree*)rl->TreeK();
 275
 276   if(trK){
 277     Int_t nlabels =0;
 278     for (Int_t i=0;i<10;i++) if (lab[i]>=0) nlabels++;
 279     if(nlabels == 0) return; // In case of no labels just exit
 280
 281
 282     Int_t ntracks = gAlice->GetMCApp()->GetNtrack();
 283
 284     for (Int_t i=0;i<10;i++){
 285       Int_t label = lab[i];
 286       if (label>=0 && label<ntracks) {
 287         TParticle *part=(TParticle*)gAlice->GetMCApp()->Particle(label);
 288
 289         if (part->P() < 0.02) {
 290           Int_t m=part->GetFirstMother();
 291           if (m<0) {
 292             continue;
 293           }
 294           if (part->GetStatusCode()>0) {
 295             continue;
 296           }
 297           lab[i]=m;
 298         }
 299         else
 300           if (part->P() < 0.12 && nlabels>3) {
 301             lab[i]=-2;
 302             nlabels--;
 303           }
 304       }
 305       else{
 306         if ( (label>ntracks||label <0) && nlabels>3) {
 307           lab[i]=-2;
 308           nlabels--;
 309         }
 310       }
 311     }
 312     if (nlabels>3){
 313       for (Int_t i=0;i<10;i++){
 314         if (nlabels>3){
 315           Int_t label = lab[i];
 316           if (label>=0 && label<ntracks) {
 317             TParticle *part=(TParticle*)gAlice->GetMCApp()->Particle(label);
 318             if (part->P() < 0.1) {
 319               lab[i]=-2;
 320               nlabels--;
 321             }
 322           }
 323         }
 324       }
 325     }
 326
 327     //compress labels -- if multi-times the same
 328     Int_t lab2[10];
 329     for (Int_t i=0;i<10;i++) lab2[i]=-2;
 330     for (Int_t i=0;i<10  ;i++){
 331       if (lab[i]<0) continue;
 332       for (Int_t j=0;j<10 &&lab2[j]!=lab[i];j++){
 333         if (lab2[j]<0) {
 334           lab2[j]= lab[i];
 335           break;
 336         }
 337       }
 338     }
 339     for (Int_t j=0;j<10;j++) lab[j]=lab2[j];
 340
 341   }
 342 }
 343
 344 //______________________________________________________________________
 345 void AliITSClusterFinder::AddLabel(Int_t lab[10], Int_t label) {
 346   //add label to the cluster
 347   AliRunLoader *rl = AliRunLoader::Instance();
 348   TTree *trK=(TTree*)rl->TreeK();
 349   if(trK){
 350     if(label<0) return; // In case of no label just exit
 351
 352     Int_t ntracks = gAlice->GetMCApp()->GetNtrack();
 353     if (label>ntracks) return;
 354     for (Int_t i=0;i<10;i++){
 355       //    if (label<0) break;
 356       if (lab[i]==label) break;
 357       if (lab[i]<0) {
 358         lab[i]= label;
 359         break;
 360       }
 361     }
 362   }
 363 }
 364
 365
 366 //______________________________________________________________________
 367 void AliITSClusterFinder::
 368 FindCluster(Int_t k,Int_t maxz,AliBin *bins,Int_t &n,Int_t *idx) {
 369   //------------------------------------------------------------
 370   // returns an array of indices of digits belonging to the cluster
 371   // (needed when the segmentation is not regular)
 372   //------------------------------------------------------------
 373   if (n<200) idx[n++]=bins[k].GetIndex();
 374   bins[k].Use();
 375
 376   if (bins[k-maxz].IsNotUsed()) FindCluster(k-maxz,maxz,bins,n,idx);
 377   if (bins[k-1   ].IsNotUsed()) FindCluster(k-1   ,maxz,bins,n,idx);
 378   if (bins[k+maxz].IsNotUsed()) FindCluster(k+maxz,maxz,bins,n,idx);
 379   if (bins[k+1   ].IsNotUsed()) FindCluster(k+1   ,maxz,bins,n,idx);
 380   /*
 381   if (bins[k-maxz-1].IsNotUsed()) FindCluster(k-maxz-1,maxz,bins,n,idx);
 382   if (bins[k-maxz+1].IsNotUsed()) FindCluster(k-maxz+1,maxz,bins,n,idx);
 383   if (bins[k+maxz-1].IsNotUsed()) FindCluster(k+maxz-1,maxz,bins,n,idx);
 384   if (bins[k+maxz+1].IsNotUsed()) FindCluster(k+maxz+1,maxz,bins,n,idx);
 385   */
 386 }
 387
 388 //______________________________________________________________________
 389 Bool_t AliITSClusterFinder::IsMaximum(Int_t k,Int_t max,const AliBin *bins) {
 390   //------------------------------------------------------------
 391   //is this a local maximum ?
 392   //------------------------------------------------------------
 393   UShort_t q=bins[k].GetQ();
 394   if (q==1023) return kFALSE;
 395   if (bins[k-max].GetQ() > q) return kFALSE;
 396   if (bins[k-1  ].GetQ() > q) return kFALSE;
 397   if (bins[k+max].GetQ() > q) return kFALSE;
 398   if (bins[k+1  ].GetQ() > q) return kFALSE;
 399   if (bins[k-max-1].GetQ() > q) return kFALSE;
 400   if (bins[k+max-1].GetQ() > q) return kFALSE;
 401   if (bins[k+max+1].GetQ() > q) return kFALSE;
 402   if (bins[k-max+1].GetQ() > q) return kFALSE;
 403   return kTRUE;
 404 }
 405
 406 //______________________________________________________________________
 407 void AliITSClusterFinder::
 408 FindPeaks(Int_t k,Int_t max,AliBin *b,Int_t *idx,UInt_t *msk,Int_t& n) {
 409   //------------------------------------------------------------
 410   //find local maxima
 411   //------------------------------------------------------------
 412   if (n<31)
 413   if (IsMaximum(k,max,b)) {
 414     idx[n]=k; msk[n]=(2<<n);
 415     n++;
 416   }
 417   b[k].SetMask(0);
 418   if (b[k-max].GetMask()&1) FindPeaks(k-max,max,b,idx,msk,n);
 419   if (b[k-1  ].GetMask()&1) FindPeaks(k-1  ,max,b,idx,msk,n);
 420   if (b[k+max].GetMask()&1) FindPeaks(k+max,max,b,idx,msk,n);
 421   if (b[k+1  ].GetMask()&1) FindPeaks(k+1  ,max,b,idx,msk,n);
 422 }
 423
 424 //______________________________________________________________________
 425 void AliITSClusterFinder::
 426 MarkPeak(Int_t k, Int_t max, AliBin *bins, UInt_t m) {
 427   //------------------------------------------------------------
 428   //mark this peak
 429   //------------------------------------------------------------
 430   UShort_t q=bins[k].GetQ();
 431
 432   bins[k].SetMask(bins[k].GetMask()|m);
 433
 434   if (bins[k-max].GetQ() <= q)
 435      if ((bins[k-max].GetMask()&m) == 0) MarkPeak(k-max,max,bins,m);
 436   if (bins[k-1  ].GetQ() <= q)
 437      if ((bins[k-1  ].GetMask()&m) == 0) MarkPeak(k-1  ,max,bins,m);
 438   if (bins[k+max].GetQ() <= q)
 439      if ((bins[k+max].GetMask()&m) == 0) MarkPeak(k+max,max,bins,m);
 440   if (bins[k+1  ].GetQ() <= q)
 441      if ((bins[k+1  ].GetMask()&m) == 0) MarkPeak(k+1  ,max,bins,m);
 442 }
 443
 444 //______________________________________________________________________
 445 void AliITSClusterFinder::
 446 MakeCluster(Int_t k,Int_t max,AliBin *bins,UInt_t m,AliITSRecPoint &c) {
 447   //------------------------------------------------------------
 448   //make cluster using digits of this peak
 449   //------------------------------------------------------------
 450   Float_t q=(Float_t)bins[k].GetQ();
 451   Int_t i=k/max, j=k-i*max;
 452
 453   c.SetQ(c.GetQ()+q);
 454   c.SetY(c.GetY()+i*q);
 455   c.SetZ(c.GetZ()+j*q);
 456   c.SetSigmaY2(c.GetSigmaY2()+i*i*q);
 457   c.SetSigmaZ2(c.GetSigmaZ2()+j*j*q);
 458
 459   bins[k].SetMask(0xFFFFFFFE);
 460
 461   if (bins[k-max].GetMask() == m) MakeCluster(k-max,max,bins,m,c);
 462   if (bins[k-1  ].GetMask() == m) MakeCluster(k-1  ,max,bins,m,c);
 463   if (bins[k+max].GetMask() == m) MakeCluster(k+max,max,bins,m,c);
 464   if (bins[k+1  ].GetMask() == m) MakeCluster(k+1  ,max,bins,m,c);
 465 }