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 <TArrayI.h>
  32 #include "AliMC.h"
  33 #include "AliLog.h"
  34
  35 ClassImp(AliITSClusterFinder)
  36
  37 extern AliRun *gAlice;
  38
  39 //----------------------------------------------------------------------
  40 AliITSClusterFinder::AliITSClusterFinder():
  41 TObject(),
  42 fModule(0),
  43 fDigits(0),
  44 fNdigits(0),
  45 fDetTypeRec(0),
  46 fClusters(0),
  47 fMap(0),
  48 fNPeaks(-1),
  49 fNModules(AliITSgeomTGeo::GetNModules()),
  50 fEvent(0),
  51 fZmin(0),
  52 fZmax(0),
  53 fXmin(0),
  54 fXmax(0),
  55 fNClusters(0),
  56 fRawID2ClusID(0)
  57 {
  58     // default cluster finder
  59     // Input:
  60     //   none.
  61     // Output:
  62     //   none.
  63     // Return:
  64     //   A default constructed AliITSCulsterFinder
  65   for(Int_t i=0; i<2200; i++){
  66     fNdet[i]=0;
  67     fNlayer[i]=0;
  68   }
  69 }
  70 //----------------------------------------------------------------------
  71 AliITSClusterFinder::AliITSClusterFinder(AliITSDetTypeRec* dettyp):
  72 TObject(),
  73 fModule(0),
  74 fDigits(0),
  75 fNdigits(0),
  76 fDetTypeRec(dettyp),
  77 fClusters(0),
  78 fMap(0),
  79 fNPeaks(-1),
  80 fNModules(AliITSgeomTGeo::GetNModules()),
  81 fEvent(0),
  82 fZmin(0),
  83 fZmax(0),
  84 fXmin(0),
  85 fXmax(0),
  86 fNClusters(0),
  87 fRawID2ClusID(0)
  88 {
  89     // default cluster finder
  90     // Standard constructor for cluster finder
  91     // Input:
  92     //   AliITSsegmentation *seg  The segmentation class to be used
  93     //   AliITSresponse     *res  The response class to be used
  94     // Output:
  95     //   none.
  96     // Return:
  97     //   A Standard constructed AliITSCulsterFinder
  98   for(Int_t i=0; i<2200; i++){
  99     fNdet[i]=0;
 100     fNlayer[i]=0;
 101   }
 102 }
 103 //----------------------------------------------------------------------
 104 AliITSClusterFinder::AliITSClusterFinder(AliITSDetTypeRec* dettyp,
 105                                          TClonesArray *digits):
 106 TObject(),
 107 fModule(0),
 108 fDigits(digits),
 109 fNdigits(0),
 110 fDetTypeRec(dettyp),
 111 fClusters(0),
 112 fMap(0),
 113 fNPeaks(-1),
 114 fNModules(AliITSgeomTGeo::GetNModules()),
 115 fEvent(0),
 116 fZmin(0),
 117 fZmax(0),
 118 fXmin(0),
 119 fXmax(0),
 120 fNClusters(0),
 121 fRawID2ClusID(0)
 122 {
 123     // default cluster finder
 124     // Standard + cluster finder constructor
 125     // Input:
 126     //   AliITSsegmentation *seg  The segmentation class to be used
 127     //   AliITSresponse     *res  The response class to be used
 128     //   TClonesArray    *digits  Array of digits to be used
 129     // Output:
 130     //   none.
 131     // Return:
 132     //   A Standard constructed AliITSCulsterFinder
 133
 134   fNdigits = fDigits->GetEntriesFast();
 135   for(Int_t i=0; i<2200; i++){
 136     fNdet[i]=0;
 137     fNlayer[i]=0;
 138   }
 139 }
 140
 141 //______________________________________________________________________
 142 AliITSClusterFinder::AliITSClusterFinder(const AliITSClusterFinder &source) :
 143   TObject(source),
 144   fModule(source.fModule),
 145   fDigits(),
 146   fNdigits(source.fNdigits),
 147   fDetTypeRec(),
 148   fClusters(),
 149   fMap(),
 150   fNPeaks(source.fNPeaks),
 151   fNModules(source.fNModules),
 152   fEvent(source.fEvent),
 153   fZmin(source.fZmin),
 154   fZmax(source.fZmax),
 155   fXmin(source.fXmin),
 156   fXmax(source.fXmax),
 157   fNClusters(source.fNClusters),
 158   fRawID2ClusID(source.fRawID2ClusID)
 159 {
 160   // Copy constructor
 161   // Copies are not allowed. The method is protected to avoid misuse.
 162   AliError("Copy constructor not allowed\n");
 163 }
 164
 165
 166 //______________________________________________________________________
 167 //AliITSClusterFinder& AliITSClusterFinder::operator=(const AliITSClusterFinder& /* source */){
 168   // Assignment operator
 169   // Assignment is not allowed. The method is protected to avoid misuse.
 170 //  Fatal("= operator","Assignment operator not allowed\n");
 171 //  return *this;
 172 //}
 173
 174 //----------------------------------------------------------------------
 175 AliITSClusterFinder::~AliITSClusterFinder(){
 176     // destructor cluster finder
 177     // Input:
 178     //   none.
 179     // Output:
 180     //   none.
 181     // Return:
 182     //   none.
 183
 184     if(fMap) {delete fMap;}
 185     // Zero local pointers. Other classes own these pointers.
 186     fMap          = 0;
 187     fDigits       = 0;
 188     fNdigits      = 0;
 189     fNPeaks       = 0;
 190     fDetTypeRec   = 0;
 191
 192 }
 193 //__________________________________________________________________________
 194 void AliITSClusterFinder::InitGeometry(){
 195  //
 196  // Initialisation of ITS geometry
 197  //
 198   Int_t mmax=AliITSgeomTGeo::GetNModules();
 199   for (Int_t m=0; m<mmax; m++) {
 200     Int_t lay,lad,det; AliITSgeomTGeo::GetModuleId(m,lay,lad,det);
 201     fNdet[m] = (lad-1)*AliITSgeomTGeo::GetNDetectors(lay) + (det-1);
 202     fNlayer[m] = lay-1;
 203   }
 204 }
 205
 206
 207
 208
 209 //______________________________________________________________________
 210 Bool_t AliITSClusterFinder::IsNeighbor(TObjArray *digs,Int_t i,Int_t n[])const{
 211     // Locagical function which checks to see if digit i has a neighbor.
 212     // If so, then it returns kTRUE and its neighbor index j.
 213     // This routine checks if the digits are side by side or one before the
 214     // other. Requires that the array of digits be in proper order.
 215     // Returns kTRUE in the following cases.
 216     //                 ji   0j   if kdiagonal  j0    0i
 217     //                 00   0i   if kdiagonal  0i    j0
 218     // Inputs:
 219     //    TObjArray *digs   Array to search for neighbors in
 220     //    Int_t      i      Index of digit for which we are searching for
 221     //                      a neighbor of.
 222     // Output:
 223     //    Int_t      j[4]   Index of one or more of the digits which is a
 224     //                      neighbor of digit a index i.
 225     // Return:
 226     //    Bool_t            kTRUE if a neighbor was found kFALSE otherwise.
 227     Int_t ix,jx,iz,jz,j;
 228     const Bool_t kdiagonal=kFALSE;
 229     Bool_t nei[4];
 230
 231     // No neighbors found if array empty.
 232     if(digs->GetEntriesFast()<=0) return kFALSE;
 233     // can not be a digit with first element or elements out or range
 234     if(i<=0 || i>=digs->GetEntriesFast()) return kFALSE;
 235
 236     for(j=0;j<4;j++){n[j] = -1;nei[j]=kFALSE;}
 237     ix = ((AliITSdigit*)(digs->At(i)))->GetCoord1();
 238     iz = ((AliITSdigit*)(digs->At(i)))->GetCoord2();
 239     for(j=0;j<i;j++){
 240         jx = ((AliITSdigit*)(digs->At(j)))->GetCoord1();
 241         jz = ((AliITSdigit*)(digs->At(j)))->GetCoord2();
 242         if(jx+1==ix && jz  ==iz){n[0] = j;nei[0] = kTRUE;}
 243         if(jx  ==ix && jz+1==iz){n[1] = j;nei[1] = kTRUE;}
 244         if(jx+1==ix && jz+1==iz){n[2] = j;nei[2] = kTRUE;}
 245         if(jx+1==ix && jz-1==iz){n[3] = j;nei[3] = kTRUE;}
 246     } // end for k
 247     if(nei[0]||nei[1]) return kTRUE;
 248     if(kdiagonal&&(nei[2]||nei[3])) return kTRUE;
 249     // no Neighbors found.
 250     return kFALSE;
 251 }
 252
 253 //______________________________________________________________________
 254 void AliITSClusterFinder::Print(ostream *os) const{
 255     //Standard output format for this class
 256     // Inputs:
 257     //    ostream *os   Output stream
 258     // Output:
 259     //    ostream *os   Output stream
 260     // Return:
 261     //    none.
 262
 263     *os << fModule<<",";
 264     *os << fNdigits<<",";
 265     *os << fNPeaks<<endl;
 266 }
 267 //______________________________________________________________________
 268 void AliITSClusterFinder::Read(istream *is)  {
 269     //Standard input for this class
 270     // Inputs:
 271     //    istream *is   Input stream
 272     // Output:
 273     //    istream *is   Input stream
 274     // Return:
 275     //    none.
 276
 277     *is >> fModule;
 278     *is >> fNdigits;
 279     *is >> fNPeaks;
 280 }
 281 //______________________________________________________________________
 282 ostream &operator<<(ostream &os,AliITSClusterFinder &source){
 283     // Standard output streaming function.
 284     // Inputs:
 285     //    ostream             *os     Output stream
 286     //    AliITSClusterFinder &source Class to be printed
 287     // Output:
 288     //    ostream             *os     Output stream
 289     // Return:
 290     //    none.
 291
 292     source.Print(&os);
 293     return os;
 294 }
 295 //______________________________________________________________________
 296 istream &operator>>(istream &is,AliITSClusterFinder &source){
 297     // Standard output streaming function.
 298     // Inputs:
 299     //    istream              *is      Input stream
 300     //     AliITSClusterFinder &source  Class to be read in.
 301     // Output:
 302     //    istream              *is      Input stream
 303     // Return:
 304     //    none.
 305
 306     source.Read(&is);
 307     return is;
 308 }
 309
 310 //______________________________________________________________________
 311 void AliITSClusterFinder::CheckLabels2(Int_t lab[10])
 312 {
 313   //------------------------------------------------------------
 314   // Tries to find mother's labels
 315   //------------------------------------------------------------
 316   AliRunLoader *rl = AliRunLoader::Instance();
 317   if(!rl) return;
 318   TTree *trK=(TTree*)rl->TreeK();
 319   if (!trK) return;
 320   //
 321   int labS[10];
 322   Int_t nlabels = 0;
 323   Int_t ntracks = gAlice->GetMCApp()->GetNtrack();
 324   for (Int_t i=0;i<10;i++) if (lab[i]>=0) labS[nlabels++] = lab[i];
 325   if (nlabels==0) return;
 326   //
 327   float mom[10];
 328   for (Int_t i=0;i<nlabels;i++) {
 329     Int_t label = labS[i];
 330     mom[i] = 0;
 331     if (label>=ntracks) continue;
 332     TParticle *part=(TParticle*)gAlice->GetMCApp()->Particle(label);
 333     mom[i] = part->P();
 334     if (part->P() < 0.02) {    // reduce soft particles from the same cluster
 335       Int_t m=part->GetFirstMother();
 336       if (m<0) continue; // primary
 337       //
 338       if (part->GetStatusCode()>0) continue;
 339       //
 340       // if the parent is within the same cluster, reassign the label to it
 341       for (int j=0;j<nlabels;j++) if (labS[j]==m) { labS[i] = m; break; }
 342     }
 343   }
 344   //
 345   if (nlabels>3) { // only 3 labels are stored in cluster, sort in decreasing momentum
 346     int ind[10],labSS[10];
 347     TMath::Sort(nlabels,mom,ind);
 348     for (int i=nlabels;i--;) labSS[i] = labS[i];
 349     for (int i=0;i<nlabels;i++) labS[i] = labSS[ind[i]];
 350   }
 351   //
 352   //compress labels -- if multi-times the same
 353   for (Int_t i=0;i<10;i++) lab[i]=-2;
 354   int nlabFin=0,j=0;
 355   for (int i=0;i<nlabels;i++) {
 356     for (j=0;j<nlabFin;j++) if (labS[i]==lab[j]) break; // the label already there
 357     if (j==nlabFin) lab[nlabFin++] = labS[i];
 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   if (fRawID2ClusID) { // RS: Register cluster id in raw words list
 489     int rwid = bins[k].GetRawID();
 490     if (fRawID2ClusID->GetSize()<=rwid) fRawID2ClusID->Set( (rwid+10)<<1 );
 491     (*fRawID2ClusID)[rwid] = fNClusters+1; // RS: store clID+1 as a reference to the cluster
 492   }
 493   if (bins[k-max].GetMask() == m) MakeCluster(k-max,max,bins,m,c);
 494   if (bins[k-1  ].GetMask() == m) MakeCluster(k-1  ,max,bins,m,c);
 495   if (bins[k+max].GetMask() == m) MakeCluster(k+max,max,bins,m,c);
 496   if (bins[k+1  ].GetMask() == m) MakeCluster(k+1  ,max,bins,m,c);
 497 }