ITS/AliITSClusterFinderV2.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-2003, 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 //            Implementation of the ITS clusterer V2 class                //
  17 //                                                                        //
  18 //          Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch            //
  19 //                                                                        //
  20 ///////////////////////////////////////////////////////////////////////////
  21
  22
  23 #include "AliRun.h"
  24 #include "AliITSClusterFinderV2.h"
  25 #include "AliITSRecPoint.h"
  26 #include "AliITSgeomTGeo.h"
  27 #include "AliITSDetTypeRec.h"
  28 #include <TParticle.h>
  29 #include "AliMC.h"
  30
  31 ClassImp(AliITSClusterFinderV2)
  32
  33 extern AliRun *gAlice;
  34
  35 AliITSClusterFinderV2::AliITSClusterFinderV2(AliITSDetTypeRec* dettyp):AliITSClusterFinder(dettyp),
  36 fNModules(AliITSgeomTGeo::GetNModules()),
  37 fEvent(0){
  38   //Default constructor
  39 }
  40 /*
  41 //______________________________________________________________________
  42 AliITSClusterFinderV2::AliITSClusterFinderV2(const AliITSClusterFinderV2 &source) : AliITSClusterFinder(source) {
  43   // Copy constructor
  44   // Copies are not allowed. The method is protected to avoid misuse.
  45   Fatal("AliITSClusterFinderV2","Copy constructor not allowed\n");
  46 }
  47 */
  48 //______________________________________________________________________
  49 //AliITSClusterFinderV2& AliITSClusterFinderV2::operator=(const AliITSClusterFinderV2& /* source */){
  50   // Assignment operator
  51   // Assignment is not allowed. The method is protected to avoid misuse.
  52   //Fatal("= operator","Assignment operator not allowed\n");
  53   //return *this;
  54 //}
  55
  56
  57 //______________________________________________________________________
  58 void AliITSClusterFinderV2::CheckLabels2(Int_t lab[10]) {
  59   //------------------------------------------------------------
  60   // Tries to find mother's labels
  61   //------------------------------------------------------------
  62   AliRunLoader *rl = AliRunLoader::GetRunLoader();
  63   TTree *trK=(TTree*)rl->TreeK();
  64
  65   if(trK){
  66     Int_t nlabels =0;
  67     for (Int_t i=0;i<10;i++) if (lab[i]>=0) nlabels++;
  68     if(nlabels == 0) return; // In case of no labels just exit
  69
  70
  71     Int_t ntracks = gAlice->GetMCApp()->GetNtrack();
  72
  73     for (Int_t i=0;i<10;i++){
  74       Int_t label = lab[i];
  75       if (label>=0 && label<ntracks) {
  76         TParticle *part=(TParticle*)gAlice->GetMCApp()->Particle(label);
  77
  78         if (part->P() < 0.02) {
  79           Int_t m=part->GetFirstMother();
  80           if (m<0) {
  81             continue;
  82           }
  83           if (part->GetStatusCode()>0) {
  84             continue;
  85           }
  86           lab[i]=m;
  87         }
  88         else
  89           if (part->P() < 0.12 && nlabels>3) {
  90             lab[i]=-2;
  91             nlabels--;
  92           }
  93       }
  94       else{
  95         if ( (label>ntracks||label <0) && nlabels>3) {
  96           lab[i]=-2;
  97           nlabels--;
  98         }
  99       }
 100     }
 101     if (nlabels>3){
 102       for (Int_t i=0;i<10;i++){
 103         if (nlabels>3){
 104           Int_t label = lab[i];
 105           if (label>=0 && label<ntracks) {
 106             TParticle *part=(TParticle*)gAlice->GetMCApp()->Particle(label);
 107             if (part->P() < 0.1) {
 108               lab[i]=-2;
 109               nlabels--;
 110             }
 111           }
 112         }
 113       }
 114     }
 115
 116     //compress labels -- if multi-times the same
 117     Int_t lab2[10];
 118     for (Int_t i=0;i<10;i++) lab2[i]=-2;
 119     for (Int_t i=0;i<10  ;i++){
 120       if (lab[i]<0) continue;
 121       for (Int_t j=0;j<10 &&lab2[j]!=lab[i];j++){
 122         if (lab2[j]<0) {
 123           lab2[j]= lab[i];
 124           break;
 125         }
 126       }
 127     }
 128     for (Int_t j=0;j<10;j++) lab[j]=lab2[j];
 129
 130   }
 131 }
 132
 133 //______________________________________________________________________
 134 void AliITSClusterFinderV2::AddLabel(Int_t lab[10], Int_t label) {
 135   //add label to the cluster
 136   AliRunLoader *rl = AliRunLoader::GetRunLoader();
 137   TTree *trK=(TTree*)rl->TreeK();
 138   if(trK){
 139     if(label<0) return; // In case of no label just exit
 140
 141     Int_t ntracks = gAlice->GetMCApp()->GetNtrack();
 142     if (label>ntracks) return;
 143     for (Int_t i=0;i<10;i++){
 144       //    if (label<0) break;
 145       if (lab[i]==label) break;
 146       if (lab[i]<0) {
 147         lab[i]= label;
 148         break;
 149       }
 150     }
 151   }
 152 }
 153
 154
 155 //______________________________________________________________________
 156 void AliITSClusterFinderV2::
 157 FindCluster(Int_t k,Int_t maxz,AliBin *bins,Int_t &n,Int_t *idx) {
 158   //------------------------------------------------------------
 159   // returns an array of indices of digits belonging to the cluster
 160   // (needed when the segmentation is not regular)
 161   //------------------------------------------------------------
 162   if (n<200) idx[n++]=bins[k].GetIndex();
 163   bins[k].Use();
 164
 165   if (bins[k-maxz].IsNotUsed()) FindCluster(k-maxz,maxz,bins,n,idx);
 166   if (bins[k-1   ].IsNotUsed()) FindCluster(k-1   ,maxz,bins,n,idx);
 167   if (bins[k+maxz].IsNotUsed()) FindCluster(k+maxz,maxz,bins,n,idx);
 168   if (bins[k+1   ].IsNotUsed()) FindCluster(k+1   ,maxz,bins,n,idx);
 169   /*
 170   if (bins[k-maxz-1].IsNotUsed()) FindCluster(k-maxz-1,maxz,bins,n,idx);
 171   if (bins[k-maxz+1].IsNotUsed()) FindCluster(k-maxz+1,maxz,bins,n,idx);
 172   if (bins[k+maxz-1].IsNotUsed()) FindCluster(k+maxz-1,maxz,bins,n,idx);
 173   if (bins[k+maxz+1].IsNotUsed()) FindCluster(k+maxz+1,maxz,bins,n,idx);
 174   */
 175 }
 176
 177 //______________________________________________________________________
 178 Bool_t AliITSClusterFinderV2::IsMaximum(Int_t k,Int_t max,const AliBin *bins) {
 179   //------------------------------------------------------------
 180   //is this a local maximum ?
 181   //------------------------------------------------------------
 182   UShort_t q=bins[k].GetQ();
 183   if (q==1023) return kFALSE;
 184   if (bins[k-max].GetQ() > q) return kFALSE;
 185   if (bins[k-1  ].GetQ() > q) return kFALSE;
 186   if (bins[k+max].GetQ() > q) return kFALSE;
 187   if (bins[k+1  ].GetQ() > q) return kFALSE;
 188   if (bins[k-max-1].GetQ() > q) return kFALSE;
 189   if (bins[k+max-1].GetQ() > q) return kFALSE;
 190   if (bins[k+max+1].GetQ() > q) return kFALSE;
 191   if (bins[k-max+1].GetQ() > q) return kFALSE;
 192   return kTRUE;
 193 }
 194
 195 //______________________________________________________________________
 196 void AliITSClusterFinderV2::
 197 FindPeaks(Int_t k,Int_t max,AliBin *b,Int_t *idx,UInt_t *msk,Int_t& n) {
 198   //------------------------------------------------------------
 199   //find local maxima
 200   //------------------------------------------------------------
 201   if (n<31)
 202   if (IsMaximum(k,max,b)) {
 203     idx[n]=k; msk[n]=(2<<n);
 204     n++;
 205   }
 206   b[k].SetMask(0);
 207   if (b[k-max].GetMask()&1) FindPeaks(k-max,max,b,idx,msk,n);
 208   if (b[k-1  ].GetMask()&1) FindPeaks(k-1  ,max,b,idx,msk,n);
 209   if (b[k+max].GetMask()&1) FindPeaks(k+max,max,b,idx,msk,n);
 210   if (b[k+1  ].GetMask()&1) FindPeaks(k+1  ,max,b,idx,msk,n);
 211 }
 212
 213 //______________________________________________________________________
 214 void AliITSClusterFinderV2::
 215 MarkPeak(Int_t k, Int_t max, AliBin *bins, UInt_t m) {
 216   //------------------------------------------------------------
 217   //mark this peak
 218   //------------------------------------------------------------
 219   UShort_t q=bins[k].GetQ();
 220
 221   bins[k].SetMask(bins[k].GetMask()|m);
 222
 223   if (bins[k-max].GetQ() <= q)
 224      if ((bins[k-max].GetMask()&m) == 0) MarkPeak(k-max,max,bins,m);
 225   if (bins[k-1  ].GetQ() <= q)
 226      if ((bins[k-1  ].GetMask()&m) == 0) MarkPeak(k-1  ,max,bins,m);
 227   if (bins[k+max].GetQ() <= q)
 228      if ((bins[k+max].GetMask()&m) == 0) MarkPeak(k+max,max,bins,m);
 229   if (bins[k+1  ].GetQ() <= q)
 230      if ((bins[k+1  ].GetMask()&m) == 0) MarkPeak(k+1  ,max,bins,m);
 231 }
 232
 233 //______________________________________________________________________
 234 void AliITSClusterFinderV2::
 235 MakeCluster(Int_t k,Int_t max,AliBin *bins,UInt_t m,AliITSRecPoint &c) {
 236   //------------------------------------------------------------
 237   //make cluster using digits of this peak
 238   //------------------------------------------------------------
 239   Float_t q=(Float_t)bins[k].GetQ();
 240   Int_t i=k/max, j=k-i*max;
 241
 242   c.SetQ(c.GetQ()+q);
 243   c.SetY(c.GetY()+i*q);
 244   c.SetZ(c.GetZ()+j*q);
 245   c.SetSigmaY2(c.GetSigmaY2()+i*i*q);
 246   c.SetSigmaZ2(c.GetSigmaZ2()+j*j*q);
 247
 248   bins[k].SetMask(0xFFFFFFFE);
 249
 250   if (bins[k-max].GetMask() == m) MakeCluster(k-max,max,bins,m,c);
 251   if (bins[k-1  ].GetMask() == m) MakeCluster(k-1  ,max,bins,m,c);
 252   if (bins[k+max].GetMask() == m) MakeCluster(k+max,max,bins,m,c);
 253   if (bins[k+1  ].GetMask() == m) MakeCluster(k+1  ,max,bins,m,c);
 254 }