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