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