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