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 "AliTPCClustersArray.h"
  29 #include "AliTPCClustersRow.h"
  30 #include "AliDigits.h"
  31 #include "AliSimDigits.h"
  32 #include "AliTPCParam.h"
  33 #include <Riostream.h>
  34 #include <TTree.h>
  35 #include "AliRunLoader.h"
  36 #include "AliLoader.h"
  37
  38 void AliTPCclusterer::FindPeaks(Int_t k,Int_t max,
  39 AliBin *b,Int_t *idx,UInt_t *msk,Int_t& n) {
  40   //find local maxima
  41   if (n<31)
  42   if (IsMaximum(k,max,b)) {
  43     idx[n]=k; msk[n]=(2<<n);
  44     n++;
  45   }
  46   b[k].SetMask(0);
  47   if (b[k-max].GetMask()&1) FindPeaks(k-max,max,b,idx,msk,n);
  48   if (b[k-1  ].GetMask()&1) FindPeaks(k-1  ,max,b,idx,msk,n);
  49   if (b[k+max].GetMask()&1) FindPeaks(k+max,max,b,idx,msk,n);
  50   if (b[k+1  ].GetMask()&1) FindPeaks(k+1  ,max,b,idx,msk,n);
  51 }
  52
  53 void AliTPCclusterer::MarkPeak(Int_t k, Int_t max, AliBin *bins, UInt_t m) {
  54   //mark this peak
  55   UShort_t q=bins[k].GetQ();
  56
  57   bins[k].SetMask(bins[k].GetMask()|m);
  58
  59   if (bins[k-max].GetQ() <= q)
  60      if ((bins[k-max].GetMask()&m) == 0) MarkPeak(k-max,max,bins,m);
  61   if (bins[k-1  ].GetQ() <= q)
  62      if ((bins[k-1  ].GetMask()&m) == 0) MarkPeak(k-1  ,max,bins,m);
  63   if (bins[k+max].GetQ() <= q)
  64      if ((bins[k+max].GetMask()&m) == 0) MarkPeak(k+max,max,bins,m);
  65   if (bins[k+1  ].GetQ() <= q)
  66      if ((bins[k+1  ].GetMask()&m) == 0) MarkPeak(k+1  ,max,bins,m);
  67 }
  68
  69 void AliTPCclusterer::MakeCluster(Int_t k,Int_t max,AliBin *bins,UInt_t m,
  70 AliTPCcluster &c) {
  71   //make cluster using digits of this peak
  72   Float_t q=(Float_t)bins[k].GetQ();
  73   Int_t i=k/max, j=k-i*max;
  74
  75   c.SetQ(c.GetQ()+q);
  76   c.SetY(c.GetY()+i*q);
  77   c.SetZ(c.GetZ()+j*q);
  78   c.SetSigmaY2(c.GetSigmaY2()+i*i*q);
  79   c.SetSigmaZ2(c.GetSigmaZ2()+j*j*q);
  80
  81   bins[k].SetMask(0xFFFFFFFE);
  82
  83   if (bins[k-max].GetMask() == m) MakeCluster(k-max,max,bins,m,c);
  84   if (bins[k-1  ].GetMask() == m) MakeCluster(k-1  ,max,bins,m,c);
  85   if (bins[k+max].GetMask() == m) MakeCluster(k+max,max,bins,m,c);
  86   if (bins[k+1  ].GetMask() == m) MakeCluster(k+1  ,max,bins,m,c);
  87 }
  88
  89 //_____________________________________________________________________________
  90 void AliTPCclusterer::Digits2Clusters(const AliTPCParam *par, AliLoader *of, Int_t eventn)
  91 {
  92   //-----------------------------------------------------------------
  93   // This is a simple cluster finder.
  94   //-----------------------------------------------------------------
  95   AliRunLoader* rl = (AliRunLoader*)of->GetEventFolder()->FindObject(AliRunLoader::fgkRunLoaderName);
  96   rl->GetEvent(eventn);
  97   TDirectory *savedir=gDirectory;
  98
  99   if (of->TreeR() == 0x0) of->MakeTree("R");
 100
 101
 102   if (of == 0x0)
 103    {
 104      cerr<<"AliTPC::Digits2Clusters(): output file not open !\n";
 105      return;
 106    }
 107
 108   const Int_t kMAXZ=par->GetMaxTBin()+2;
 109
 110
 111   TTree *t = (TTree *)of->TreeD();
 112
 113   if (!t) {
 114     cerr<<"Input tree with digits not found"<<endl;
 115     return;
 116   }
 117
 118   AliSimDigits digarr, *dummy=&digarr;
 119   t->GetBranch("Segment")->SetAddress(&dummy);
 120   Stat_t nentries = t->GetEntries();
 121
 122   cout<<"Got "<<nentries<<" from TreeD"<<endl;
 123
 124 //  ((AliTPCParam*)par)->Write(par->GetTitle());
 125
 126   AliTPCClustersArray carray;
 127   carray.Setup(par);
 128   carray.SetClusterType("AliTPCcluster");
 129
 130   TTree* treeR = of->TreeR();
 131   carray.MakeTree(treeR);
 132
 133
 134
 135   Int_t nclusters=0;
 136
 137   for (Int_t n=0; n<nentries; n++)
 138    {
 139
 140     Int_t sec, row;
 141     t->GetEvent(n);
 142
 143     if (!par->AdjustSectorRow(digarr.GetID(),sec,row)) {
 144        cerr<<"AliTPC warning: invalid segment ID ! "<<digarr.GetID()<<endl;
 145        continue;
 146     }
 147
 148     AliTPCClustersRow *clrow=carray.CreateRow(sec,row);
 149
 150     Float_t rx=par->GetPadRowRadii(sec,row);
 151
 152     Int_t npads, sign;
 153     {
 154        const Int_t kNIS=par->GetNInnerSector(), kNOS=par->GetNOuterSector();
 155        if (sec < kNIS) {
 156           npads = par->GetNPadsLow(row);
 157           sign = (sec < kNIS/2) ? 1 : -1;
 158        } else {
 159           npads = par->GetNPadsUp(row);
 160           sign = ((sec-kNIS) < kNOS/2) ? 1 : -1;
 161        }
 162     }
 163
 164     const Int_t kMAXBIN=kMAXZ*(npads+2);
 165     AliBin *bins=new AliBin[kMAXBIN];
 166     for (Int_t ii=0;ii<kMAXBIN;ii++) {
 167        bins[ii].SetQ(0); bins[ii].SetMask(0xFFFFFFFE);
 168     }
 169
 170     digarr.First();
 171     do {
 172        Short_t dig=digarr.CurrentDigit();
 173        if (dig<=par->GetZeroSup()) continue;
 174        Int_t j=digarr.CurrentRow()+1, i=digarr.CurrentColumn()+1;
 175        bins[i*kMAXZ+j].SetQ(dig);
 176        bins[i*kMAXZ+j].SetMask(1);
 177     } while (digarr.Next());
 178
 179     Int_t ncl=0;
 180     for (Int_t i=0; i<kMAXBIN; i++) {
 181       if ((bins[i].GetMask()&1) == 0) continue;
 182       Int_t idx[32]; UInt_t msk[32]; Int_t npeaks=0;
 183       FindPeaks(i, kMAXZ, bins, idx, msk, npeaks);
 184
 185       if (npeaks>30) continue;
 186
 187       Int_t k,l;
 188       for (k=0; k<npeaks-1; k++){//mark adjacent peaks
 189         if (idx[k] < 0) continue; //this peak is already removed
 190         for (l=k+1; l<npeaks; l++) {
 191            if (idx[l] < 0) continue; //this peak is already removed
 192            Int_t ki=idx[k]/kMAXZ, kj=idx[k] - ki*kMAXZ;
 193            Int_t li=idx[l]/kMAXZ, lj=idx[l] - li*kMAXZ;
 194            Int_t di=TMath::Abs(ki - li);
 195            Int_t dj=TMath::Abs(kj - lj);
 196            if (di>1 || dj>1) continue;
 197            if (bins[idx[k]].GetQ() > bins[idx[l]].GetQ()) {
 198               msk[l]=msk[k];
 199               idx[l]*=-1;
 200            } else {
 201               msk[k]=msk[l];
 202               idx[k]*=-1;
 203               break;
 204            }
 205         }
 206       }
 207
 208       for (k=0; k<npeaks; k++) {
 209         MarkPeak(TMath::Abs(idx[k]), kMAXZ, bins, msk[k]);
 210       }
 211
 212       for (k=0; k<npeaks; k++) {
 213          if (idx[k] < 0) continue; //removed peak
 214          AliTPCcluster c;
 215          MakeCluster(idx[k], kMAXZ, bins, msk[k], c);
 216          if (c.GetQ() < 5) continue; //noise cluster
 217          c.SetY(c.GetY()/c.GetQ());
 218          c.SetZ(c.GetZ()/c.GetQ());
 219
 220          Float_t s2 = c.GetSigmaY2()/c.GetQ() - c.GetY()*c.GetY();
 221          Float_t w=par->GetPadPitchWidth(sec);
 222          c.SetSigmaY2((s2 + 1./12.)*w*w);
 223          if (s2 != 0.) {
 224            c.SetSigmaY2(c.GetSigmaY2()*0.108);
 225            if (sec<par->GetNInnerSector()) c.SetSigmaY2(c.GetSigmaY2()*2.07);
 226          }
 227
 228          s2 = c.GetSigmaZ2()/c.GetQ() - c.GetZ()*c.GetZ();
 229          w=par->GetZWidth();
 230          c.SetSigmaZ2((s2 + 1./12.)*w*w);
 231          if (s2 != 0.) {
 232            c.SetSigmaZ2(c.GetSigmaZ2()*0.169);
 233            if (sec<par->GetNInnerSector()) c.SetSigmaZ2(c.GetSigmaZ2()*1.77);
 234          }
 235
 236          c.SetY((c.GetY() - 0.5 - 0.5*npads)*par->GetPadPitchWidth(sec));
 237          c.SetZ(par->GetZWidth()*(c.GetZ()-1));
 238          c.SetZ(c.GetZ() - 3.*par->GetZSigma()); // PASA delay
 239          c.SetZ(sign*(par->GetZLength() - c.GetZ()));
 240
 241          if (rx<230./250.*TMath::Abs(c.GetZ())) continue;
 242
 243          Int_t ki=idx[k]/kMAXZ, kj=idx[k] - ki*kMAXZ;
 244          c.SetLabel(digarr.GetTrackID(kj-1,ki-1,0),0);
 245          c.SetLabel(digarr.GetTrackID(kj-1,ki-1,1),1);
 246          c.SetLabel(digarr.GetTrackID(kj-1,ki-1,2),2);
 247
 248          c.SetQ(bins[idx[k]].GetQ());
 249
 250          if (ki==1 || ki==npads || kj==1 || kj==kMAXZ-2) {
 251            c.SetSigmaY2(c.GetSigmaY2()*25.);
 252            c.SetSigmaZ2(c.GetSigmaZ2()*4.);
 253          }
 254          clrow->InsertCluster(&c); ncl++;
 255       }
 256     }
 257     carray.StoreRow(sec,row);
 258     carray.ClearRow(sec,row);
 259
 260     //cerr<<"sector, row, compressed digits, clusters: "
 261     //<<sec<<' '<<row<<' '<<digarr.GetSize()<<' '<<ncl<<"                  \r";
 262
 263     nclusters+=ncl;
 264
 265     delete[] bins;
 266   }
 267
 268   cerr<<"Number of found clusters : "<<nclusters<<"                        \n";
 269
 270   of->WriteRecPoints("OVERWRITE");
 271
 272   savedir->cd();
 273
 274 //  delete t;  //Thanks to Mariana Bondila
 275 }
 276