TPC/AliTPCclusterer.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 /* $Id$ */
  17
  18 //-------------------------------------------------------
  19 //          Implementation of the TPC clusterer
  20 //
  21 //   Origin: Jouri Belikov, CERN, Jouri.Belikov@cern.ch
  22 //-------------------------------------------------------
  23
  24 #include "AliTPCclusterer.h"
  25 #include "AliTPCcluster.h"
  26 #include <TObjArray.h>
  27 #include <TFile.h>
  28 #include "AliDigits.h"
  29 #include "AliSimDigits.h"
  30 #include "AliTPCParam.h"
  31 #include "AliTPCClustersRow.h"
  32 #include <TTree.h>
  33
  34 ClassImp(AliTPCclusterer)
  35
  36
  37 //____________________________________________________
  38 AliTPCclusterer::AliTPCclusterer(const AliTPCclusterer &param)
  39                 :TObject(),
  40                  fPar(0)
  41 {
  42   //
  43   // dummy
  44   //
  45   fPar = param.fPar;
  46 }
  47 //-----------------------------------------------------
  48 AliTPCclusterer & AliTPCclusterer::operator =(const AliTPCclusterer & param)
  49 {
  50   //
  51   // assignment operator - dummy
  52   //
  53   if (this == &param) return (*this);
  54
  55   fPar = param.fPar;
  56   return (*this);
  57 }
  58 //____________________________________________________
  59 void AliTPCclusterer::FindPeaks(Int_t k,Int_t max,
  60 AliBin *b,Int_t *idx,UInt_t *msk,Int_t& n) {
  61   //find local maxima
  62   if (n<31)
  63   if (IsMaximum(k,max,b)) {
  64     idx[n]=k; msk[n]=(2<<n);
  65     n++;
  66   }
  67   b[k].SetMask(0);
  68   if (b[k-max].GetMask()&1) FindPeaks(k-max,max,b,idx,msk,n);
  69   if (b[k-1  ].GetMask()&1) FindPeaks(k-1  ,max,b,idx,msk,n);
  70   if (b[k+max].GetMask()&1) FindPeaks(k+max,max,b,idx,msk,n);
  71   if (b[k+1  ].GetMask()&1) FindPeaks(k+1  ,max,b,idx,msk,n);
  72 }
  73
  74 void AliTPCclusterer::MarkPeak(Int_t k, Int_t max, AliBin *bins, UInt_t m) {
  75   //mark this peak
  76   UShort_t q=bins[k].GetQ();
  77
  78   bins[k].SetMask(bins[k].GetMask()|m);
  79
  80   if (bins[k-max].GetQ() <= q)
  81      if ((bins[k-max].GetMask()&m) == 0) MarkPeak(k-max,max,bins,m);
  82   if (bins[k-1  ].GetQ() <= q)
  83      if ((bins[k-1  ].GetMask()&m) == 0) MarkPeak(k-1  ,max,bins,m);
  84   if (bins[k+max].GetQ() <= q)
  85      if ((bins[k+max].GetMask()&m) == 0) MarkPeak(k+max,max,bins,m);
  86   if (bins[k+1  ].GetQ() <= q)
  87      if ((bins[k+1  ].GetMask()&m) == 0) MarkPeak(k+1  ,max,bins,m);
  88 }
  89
  90 void AliTPCclusterer::MakeCluster(Int_t k,Int_t max,AliBin *bins,UInt_t m,
  91 AliTPCcluster &c) {
  92   //make cluster using digits of this peak
  93   Float_t q=(Float_t)bins[k].GetQ();
  94   Int_t i=k/max, j=k-i*max;
  95
  96   c.SetQ(c.GetQ()+q);
  97   c.SetY(c.GetY()+i*q);
  98   c.SetZ(c.GetZ()+j*q);
  99   c.SetSigmaY2(c.GetSigmaY2()+i*i*q);
 100   c.SetSigmaZ2(c.GetSigmaZ2()+j*j*q);
 101
 102   bins[k].SetMask(0xFFFFFFFE);
 103
 104   if (bins[k-max].GetMask() == m) MakeCluster(k-max,max,bins,m,c);
 105   if (bins[k-1  ].GetMask() == m) MakeCluster(k-1  ,max,bins,m,c);
 106   if (bins[k+max].GetMask() == m) MakeCluster(k+max,max,bins,m,c);
 107   if (bins[k+1  ].GetMask() == m) MakeCluster(k+1  ,max,bins,m,c);
 108 }
 109
 110 //_____________________________________________________________________________
 111 Int_t AliTPCclusterer::Digits2Clusters(TTree *dTree, TTree *cTree) {
 112   //-----------------------------------------------------------------
 113   // This is a simple cluster finder.
 114   //-----------------------------------------------------------------
 115   TBranch *branch=dTree->GetBranch("Segment");
 116   if (!branch) {
 117      Error("Digits2Cluster","Can't get the branch !");
 118      return 1;
 119   }
 120   AliSimDigits digarr, *dummy=&digarr;
 121   branch->SetAddress(&dummy);
 122
 123   AliTPCClustersRow ddd,*clrow=&ddd;
 124   clrow = new AliTPCClustersRow("AliTPCcluster");
 125   cTree->Branch("Segment","AliTPCClustersRow",&clrow,32000,200);
 126
 127   const Int_t kMAXZ=fPar->GetMaxTBin()+2;
 128
 129   Int_t nclusters=0;
 130
 131   Int_t nentries = (Int_t)dTree->GetEntries();
 132   for (Int_t n=0; n<nentries; n++) {
 133
 134     Int_t sec, row;
 135     dTree->GetEvent(n);
 136
 137     if (!fPar->AdjustSectorRow(digarr.GetID(),sec,row)) {
 138        Error("Digits2Clusters","!nvalid segment ID ! %d",digarr.GetID());
 139        continue;
 140     }
 141
 142     clrow->SetID(digarr.GetID());
 143
 144     cTree->GetBranch("Segment")->SetAddress(&clrow);
 145
 146     Float_t rx=fPar->GetPadRowRadii(sec,row);
 147
 148     Int_t npads, sign;
 149     {
 150        const Int_t kNIS=fPar->GetNInnerSector(), kNOS=fPar->GetNOuterSector();
 151        if (sec < kNIS) {
 152           npads = fPar->GetNPadsLow(row);
 153           sign = (sec < kNIS/2) ? 1 : -1;
 154        } else {
 155           npads = fPar->GetNPadsUp(row);
 156           sign = ((sec-kNIS) < kNOS/2) ? 1 : -1;
 157        }
 158     }
 159
 160     const Int_t kMAXBIN=kMAXZ*(npads+2);
 161     AliBin *bins=new AliBin[kMAXBIN];
 162     for (Int_t ii=0;ii<kMAXBIN;ii++) {
 163        bins[ii].SetQ(0); bins[ii].SetMask(0xFFFFFFFE);
 164     }
 165
 166     digarr.First();
 167     do {
 168        Short_t dig=digarr.CurrentDigit();
 169        if (dig<=fPar->GetZeroSup()) continue;
 170        Int_t j=digarr.CurrentRow()+1, i=digarr.CurrentColumn()+1;
 171        bins[i*kMAXZ+j].SetQ(dig);
 172        bins[i*kMAXZ+j].SetMask(1);
 173     } while (digarr.Next());
 174
 175     Int_t ncl=0;
 176     for (Int_t i=0; i<kMAXBIN; i++) {
 177       if ((bins[i].GetMask()&1) == 0) continue;
 178       Int_t idx[32]; UInt_t msk[32]; Int_t npeaks=0;
 179       FindPeaks(i, kMAXZ, bins, idx, msk, npeaks);
 180
 181       if (npeaks>30) continue;
 182
 183       Int_t k,l;
 184       for (k=0; k<npeaks-1; k++){//mark adjacent peaks
 185         if (idx[k] < 0) continue; //this peak is already removed
 186         for (l=k+1; l<npeaks; l++) {
 187            if (idx[l] < 0) continue; //this peak is already removed
 188            Int_t ki=idx[k]/kMAXZ, kj=idx[k] - ki*kMAXZ;
 189            Int_t li=idx[l]/kMAXZ, lj=idx[l] - li*kMAXZ;
 190            Int_t di=TMath::Abs(ki - li);
 191            Int_t dj=TMath::Abs(kj - lj);
 192            if (di>1 || dj>1) continue;
 193            if (bins[idx[k]].GetQ() > bins[idx[l]].GetQ()) {
 194               msk[l]=msk[k];
 195               idx[l]*=-1;
 196            } else {
 197               msk[k]=msk[l];
 198               idx[k]*=-1;
 199               break;
 200            }
 201         }
 202       }
 203
 204       for (k=0; k<npeaks; k++) {
 205         MarkPeak(TMath::Abs(idx[k]), kMAXZ, bins, msk[k]);
 206       }
 207
 208       for (k=0; k<npeaks; k++) {
 209          if (idx[k] < 0) continue; //removed peak
 210          AliTPCcluster c;
 211          MakeCluster(idx[k], kMAXZ, bins, msk[k], c);
 212          if (c.GetQ() < 5) continue; //noise cluster
 213          c.SetY(c.GetY()/c.GetQ());
 214          c.SetZ(c.GetZ()/c.GetQ());
 215
 216          Float_t s2 = c.GetSigmaY2()/c.GetQ() - c.GetY()*c.GetY();
 217          Float_t w=fPar->GetPadPitchWidth(sec);
 218          c.SetSigmaY2((s2 + 1./12.)*w*w);
 219          if (s2 != 0.) {
 220            c.SetSigmaY2(c.GetSigmaY2()*0.108);
 221            if (sec<fPar->GetNInnerSector()) c.SetSigmaY2(c.GetSigmaY2()*2.07);
 222          }
 223
 224          s2 = c.GetSigmaZ2()/c.GetQ() - c.GetZ()*c.GetZ();
 225          w=fPar->GetZWidth();
 226          c.SetSigmaZ2((s2 + 1./12.)*w*w);
 227          if (s2 != 0.) {
 228            c.SetSigmaZ2(c.GetSigmaZ2()*0.169);
 229            if (sec<fPar->GetNInnerSector()) c.SetSigmaZ2(c.GetSigmaZ2()*1.77);
 230          }
 231
 232          c.SetY((c.GetY() - 0.5 - 0.5*npads)*fPar->GetPadPitchWidth(sec));
 233          c.SetZ(fPar->GetZWidth()*(c.GetZ()-1));
 234          c.SetZ(c.GetZ() - 3.*fPar->GetZSigma()); // PASA delay
 235          c.SetZ(sign*(fPar->GetZLength(sec) - c.GetZ()));
 236
 237          if (rx<230./250.*TMath::Abs(c.GetZ())) continue;
 238
 239          Int_t ki=idx[k]/kMAXZ, kj=idx[k] - ki*kMAXZ;
 240          c.SetLabel(digarr.GetTrackID(kj-1,ki-1,0),0);
 241          c.SetLabel(digarr.GetTrackID(kj-1,ki-1,1),1);
 242          c.SetLabel(digarr.GetTrackID(kj-1,ki-1,2),2);
 243
 244          c.SetQ(bins[idx[k]].GetQ());
 245
 246          if (ki==1 || ki==npads || kj==1 || kj==kMAXZ-2) {
 247            c.SetSigmaY2(c.GetSigmaY2()*25.);
 248            c.SetSigmaZ2(c.GetSigmaZ2()*4.);
 249          }
 250          clrow->InsertCluster(&c); ncl++;
 251       }
 252     }
 253     cTree->Fill();
 254
 255     clrow->GetArray()->Clear();
 256
 257     nclusters+=ncl;
 258
 259     delete[] bins;
 260   }
 261
 262   Info("Digits2Cluster","Number of found clusters : %d",nclusters);
 263
 264   return 0;
 265 }
 266