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