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 "AliITSclusterV2.h"
  27 #include "AliITSsimulationFastPoints.h"
  28 #include "AliITS.h"
  29 #include "AliITSgeom.h"
  30 #include <TParticle.h>
  31 #include "AliMC.h"
  32
  33 ClassImp(AliITSClusterFinderV2)
  34
  35 extern AliRun *gAlice;
  36
  37 AliITSClusterFinderV2::AliITSClusterFinderV2():AliITSClusterFinder(){
  38
  39   //Default constructor
  40   fEvent = 0;
  41   fModule = 0;
  42
  43   AliITSgeom *geom=(AliITSgeom*)fITS->GetITSgeom();
  44
  45   fNModules = geom->GetIndexMax();
  46 }
  47
  48 //______________________________________________________________________
  49 AliITSClusterFinderV2::AliITSClusterFinderV2(const AliITSClusterFinderV2 &source) : AliITSClusterFinder(source) {
  50   // Copy constructor
  51   // Copies are not allowed. The method is protected to avoid misuse.
  52   Fatal("AliITSClusterFinderV2","Copy constructor not allowed\n");
  53 }
  54
  55 //______________________________________________________________________
  56 AliITSClusterFinderV2& AliITSClusterFinderV2::operator=(const AliITSClusterFinderV2& /* source */){
  57   // Assignment operator
  58   // Assignment is not allowed. The method is protected to avoid misuse.
  59   Fatal("= operator","Assignment operator not allowed\n");
  60   return *this;
  61 }
  62
  63
  64 //______________________________________________________________________
  65 void AliITSClusterFinderV2::CheckLabels2(Int_t lab[10]) {
  66   //------------------------------------------------------------
  67   // Tries to find mother's labels
  68   //------------------------------------------------------------
  69   Int_t nlabels =0;
  70   for (Int_t i=0;i<10;i++) if (lab[i]>=0) nlabels++;
  71   if(nlabels == 0) return; // In case of no labels just exit
  72
  73
  74   Int_t ntracks = gAlice->GetMCApp()->GetNtrack();
  75
  76   for (Int_t i=0;i<10;i++){
  77     Int_t label = lab[i];
  78     if (label>=0 && label<ntracks) {
  79       TParticle *part=(TParticle*)gAlice->GetMCApp()->Particle(label);
  80
  81       if (part->P() < 0.02) {
  82           Int_t m=part->GetFirstMother();
  83           if (m<0) {
  84             continue;
  85           }
  86           if (part->GetStatusCode()>0) {
  87             continue;
  88           }
  89           lab[i]=m;
  90       }
  91       else
  92         if (part->P() < 0.12 && nlabels>3) {
  93           lab[i]=-2;
  94           nlabels--;
  95         }
  96     }
  97     else{
  98       if ( (label>ntracks||label <0) && nlabels>3) {
  99         lab[i]=-2;
 100         nlabels--;
 101       }
 102     }
 103   }
 104   if (nlabels>3){
 105     for (Int_t i=0;i<10;i++){
 106       if (nlabels>3){
 107         Int_t label = lab[i];
 108         if (label>=0 && label<ntracks) {
 109           TParticle *part=(TParticle*)gAlice->GetMCApp()->Particle(label);
 110           if (part->P() < 0.1) {
 111             lab[i]=-2;
 112             nlabels--;
 113           }
 114         }
 115       }
 116     }
 117   }
 118
 119   //compress labels -- if multi-times the same
 120   Int_t lab2[10];
 121   for (Int_t i=0;i<10;i++) lab2[i]=-2;
 122   for (Int_t i=0;i<10  ;i++){
 123     if (lab[i]<0) continue;
 124     for (Int_t j=0;j<10 &&lab2[j]!=lab[i];j++){
 125       if (lab2[j]<0) {
 126         lab2[j]= lab[i];
 127         break;
 128       }
 129     }
 130   }
 131   for (Int_t j=0;j<10;j++) lab[j]=lab2[j];
 132
 133 }
 134
 135 //______________________________________________________________________
 136 void AliITSClusterFinderV2::AddLabel(Int_t lab[10], Int_t label) {
 137
 138   //add label to the cluster
 139
 140   if(label<0) return; // In case of no label just exit
 141
 142   Int_t ntracks = gAlice->GetMCApp()->GetNtrack();
 143   if (label>ntracks) return;
 144   for (Int_t i=0;i<10;i++){
 145     //    if (label<0) break;
 146     if (lab[i]==label) break;
 147     if (lab[i]<0) {
 148       lab[i]= label;
 149       break;
 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,AliITSclusterV2 &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 }
 255
 256 //______________________________________________________________________
 257 Int_t AliITSClusterFinderV2::Hits2Clusters(TTree *hTree, TTree *cTree) {
 258   //------------------------------------------------------------
 259   // This function creates ITS clusters
 260   //------------------------------------------------------------
 261
 262   AliITSgeom *geom=fITS->GetITSgeom();
 263   Int_t mmax=geom->GetIndexMax();
 264
 265   fITS->InitModules(-1,mmax);
 266   fITS->FillModules(hTree,0);
 267
 268   TClonesArray *clusters=new TClonesArray("AliITSclusterV2",1000);
 269   TBranch *branch=cTree->GetBranch("Clusters");
 270   if (!branch) cTree->Branch("Clusters",&clusters);
 271   else branch->SetAddress(&clusters);
 272
 273   static TClonesArray *points=fITS->RecPoints();
 274   AliITSsimulationFastPoints sim;
 275   Int_t ncl=0;
 276   for (Int_t m=0; m<mmax; m++) {
 277     AliITSmodule *mod=fITS->GetModule(m);
 278     sim.CreateFastRecPoints(mod,m,gRandom);
 279
 280     RecPoints2Clusters(points, m, clusters);
 281     fITS->ResetRecPoints();
 282
 283     ncl+=clusters->GetEntriesFast();
 284     cTree->Fill();
 285     clusters->Clear();
 286   }
 287
 288   Info("Hits2Clusters","Number of found fast clusters : %d",ncl);
 289
 290   //cTree->Write();
 291
 292   delete clusters;
 293
 294   return 0;
 295 }
 296