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