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   TTree *trK=(TTree*)rl->TreeK();
 274
 275   if(trK){
 276     Int_t nlabels =0;
 277     for (Int_t i=0;i<10;i++) if (lab[i]>=0) nlabels++;
 278     if(nlabels == 0) return; // In case of no labels just exit
 279
 280
 281     Int_t ntracks = gAlice->GetMCApp()->GetNtrack();
 282
 283     for (Int_t i=0;i<10;i++){
 284       Int_t label = lab[i];
 285       if (label>=0 && label<ntracks) {
 286         TParticle *part=(TParticle*)gAlice->GetMCApp()->Particle(label);
 287
 288         if (part->P() < 0.02) {
 289           Int_t m=part->GetFirstMother();
 290           if (m<0) {
 291             continue;
 292           }
 293           if (part->GetStatusCode()>0) {
 294             continue;
 295           }
 296           lab[i]=m;
 297         }
 298         else
 299           if (part->P() < 0.12 && nlabels>3) {
 300             lab[i]=-2;
 301             nlabels--;
 302           }
 303       }
 304       else{
 305         if ( (label>ntracks||label <0) && nlabels>3) {
 306           lab[i]=-2;
 307           nlabels--;
 308         }
 309       }
 310     }
 311     if (nlabels>3){
 312       for (Int_t i=0;i<10;i++){
 313         if (nlabels>3){
 314           Int_t label = lab[i];
 315           if (label>=0 && label<ntracks) {
 316             TParticle *part=(TParticle*)gAlice->GetMCApp()->Particle(label);
 317             if (part->P() < 0.1) {
 318               lab[i]=-2;
 319               nlabels--;
 320             }
 321           }
 322         }
 323       }
 324     }
 325
 326     //compress labels -- if multi-times the same
 327     Int_t lab2[10];
 328     for (Int_t i=0;i<10;i++) lab2[i]=-2;
 329     for (Int_t i=0;i<10  ;i++){
 330       if (lab[i]<0) continue;
 331       for (Int_t j=0;j<10 &&lab2[j]!=lab[i];j++){
 332         if (lab2[j]<0) {
 333           lab2[j]= lab[i];
 334           break;
 335         }
 336       }
 337     }
 338     for (Int_t j=0;j<10;j++) lab[j]=lab2[j];
 339
 340   }
 341 }
 342
 343 //______________________________________________________________________
 344 void AliITSClusterFinder::AddLabel(Int_t lab[10], Int_t label) {
 345   //add label to the cluster
 346   AliRunLoader *rl = AliRunLoader::Instance();
 347   TTree *trK=(TTree*)rl->TreeK();
 348   if(trK){
 349     if(label<0) return; // In case of no label just exit
 350
 351     Int_t ntracks = gAlice->GetMCApp()->GetNtrack();
 352     if (label>ntracks) return;
 353     for (Int_t i=0;i<10;i++){
 354       //    if (label<0) break;
 355       if (lab[i]==label) break;
 356       if (lab[i]<0) {
 357         lab[i]= label;
 358         break;
 359       }
 360     }
 361   }
 362 }
 363
 364
 365 //______________________________________________________________________
 366 void AliITSClusterFinder::
 367 FindCluster(Int_t k,Int_t maxz,AliBin *bins,Int_t &n,Int_t *idx) {
 368   //------------------------------------------------------------
 369   // returns an array of indices of digits belonging to the cluster
 370   // (needed when the segmentation is not regular)
 371   //------------------------------------------------------------
 372   if (n<200) idx[n++]=bins[k].GetIndex();
 373   bins[k].Use();
 374
 375   if (bins[k-maxz].IsNotUsed()) FindCluster(k-maxz,maxz,bins,n,idx);
 376   if (bins[k-1   ].IsNotUsed()) FindCluster(k-1   ,maxz,bins,n,idx);
 377   if (bins[k+maxz].IsNotUsed()) FindCluster(k+maxz,maxz,bins,n,idx);
 378   if (bins[k+1   ].IsNotUsed()) FindCluster(k+1   ,maxz,bins,n,idx);
 379   /*
 380   if (bins[k-maxz-1].IsNotUsed()) FindCluster(k-maxz-1,maxz,bins,n,idx);
 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   */
 385 }
 386
 387 //______________________________________________________________________
 388 Bool_t AliITSClusterFinder::IsMaximum(Int_t k,Int_t max,const AliBin *bins) {
 389   //------------------------------------------------------------
 390   //is this a local maximum ?
 391   //------------------------------------------------------------
 392   UShort_t q=bins[k].GetQ();
 393   if (q==1023) return kFALSE;
 394   if (bins[k-max].GetQ() > q) return kFALSE;
 395   if (bins[k-1  ].GetQ() > q) return kFALSE;
 396   if (bins[k+max].GetQ() > q) return kFALSE;
 397   if (bins[k+1  ].GetQ() > q) return kFALSE;
 398   if (bins[k-max-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   return kTRUE;
 403 }
 404
 405 //______________________________________________________________________
 406 void AliITSClusterFinder::
 407 FindPeaks(Int_t k,Int_t max,AliBin *b,Int_t *idx,UInt_t *msk,Int_t& n) {
 408   //------------------------------------------------------------
 409   //find local maxima
 410   //------------------------------------------------------------
 411   if (n<31)
 412   if (IsMaximum(k,max,b)) {
 413     idx[n]=k; msk[n]=(2<<n);
 414     n++;
 415   }
 416   b[k].SetMask(0);
 417   if (b[k-max].GetMask()&1) FindPeaks(k-max,max,b,idx,msk,n);
 418   if (b[k-1  ].GetMask()&1) FindPeaks(k-1  ,max,b,idx,msk,n);
 419   if (b[k+max].GetMask()&1) FindPeaks(k+max,max,b,idx,msk,n);
 420   if (b[k+1  ].GetMask()&1) FindPeaks(k+1  ,max,b,idx,msk,n);
 421 }
 422
 423 //______________________________________________________________________
 424 void AliITSClusterFinder::
 425 MarkPeak(Int_t k, Int_t max, AliBin *bins, UInt_t m) {
 426   //------------------------------------------------------------
 427   //mark this peak
 428   //------------------------------------------------------------
 429   UShort_t q=bins[k].GetQ();
 430
 431   bins[k].SetMask(bins[k].GetMask()|m);
 432
 433   if (bins[k-max].GetQ() <= q)
 434      if ((bins[k-max].GetMask()&m) == 0) MarkPeak(k-max,max,bins,m);
 435   if (bins[k-1  ].GetQ() <= q)
 436      if ((bins[k-1  ].GetMask()&m) == 0) MarkPeak(k-1  ,max,bins,m);
 437   if (bins[k+max].GetQ() <= q)
 438      if ((bins[k+max].GetMask()&m) == 0) MarkPeak(k+max,max,bins,m);
 439   if (bins[k+1  ].GetQ() <= q)
 440      if ((bins[k+1  ].GetMask()&m) == 0) MarkPeak(k+1  ,max,bins,m);
 441 }
 442
 443 //______________________________________________________________________
 444 void AliITSClusterFinder::
 445 MakeCluster(Int_t k,Int_t max,AliBin *bins,UInt_t m,AliITSRecPoint &c) {
 446   //------------------------------------------------------------
 447   //make cluster using digits of this peak
 448   //------------------------------------------------------------
 449   Float_t q=(Float_t)bins[k].GetQ();
 450   Int_t i=k/max, j=k-i*max;
 451
 452   c.SetQ(c.GetQ()+q);
 453   c.SetY(c.GetY()+i*q);
 454   c.SetZ(c.GetZ()+j*q);
 455   c.SetSigmaY2(c.GetSigmaY2()+i*i*q);
 456   c.SetSigmaZ2(c.GetSigmaZ2()+j*j*q);
 457
 458   bins[k].SetMask(0xFFFFFFFE);
 459
 460   if (bins[k-max].GetMask() == m) MakeCluster(k-max,max,bins,m,c);
 461   if (bins[k-1  ].GetMask() == m) MakeCluster(k-1  ,max,bins,m,c);
 462   if (bins[k+max].GetMask() == m) MakeCluster(k+max,max,bins,m,c);
 463   if (bins[k+1  ].GetMask() == m) MakeCluster(k+1  ,max,bins,m,c);
 464 }