Removing obsolete macros which used the AliTPCtracker
[u/mrichter/AliRoot.git] / TPC / AliTPCclusterer.cxx
CommitLineData
73042f01 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
88cb7938 16/* $Id$ */
73042f01 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>
73042f01 28#include "AliDigits.h"
29#include "AliSimDigits.h"
30#include "AliTPCParam.h"
d7cbafb5 31#include "AliTPCClustersRow.h"
c3c63118 32#include <TTree.h>
c630aafd 33
34ClassImp(AliTPCclusterer)
35
73042f01 36
179c6296 37//____________________________________________________
38AliTPCclusterer::AliTPCclusterer(const AliTPCclusterer &param)
39 :TObject(),
40 fPar(0)
41{
42 //
43 // dummy
44 //
45 fPar = param.fPar;
46}
47//-----------------------------------------------------
48AliTPCclusterer & AliTPCclusterer::operator =(const AliTPCclusterer & param)
49{
50 //
51 // assignment operator - dummy
52 //
04420071 53 if (this == &param) return (*this);
54
55 fPar = param.fPar;
179c6296 56 return (*this);
57}
58//____________________________________________________
73042f01 59void AliTPCclusterer::FindPeaks(Int_t k,Int_t max,
60AliBin *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
74void 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
90void AliTPCclusterer::MakeCluster(Int_t k,Int_t max,AliBin *bins,UInt_t m,
91AliTPCcluster &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//_____________________________________________________________________________
c630aafd 111Int_t AliTPCclusterer::Digits2Clusters(TTree *dTree, TTree *cTree) {
73042f01 112 //-----------------------------------------------------------------
113 // This is a simple cluster finder.
114 //-----------------------------------------------------------------
c630aafd 115 TBranch *branch=dTree->GetBranch("Segment");
116 if (!branch) {
117 Error("Digits2Cluster","Can't get the branch !");
118 return 1;
6e4d905e 119 }
73042f01 120 AliSimDigits digarr, *dummy=&digarr;
c630aafd 121 branch->SetAddress(&dummy);
88cb7938 122
d7cbafb5 123 AliTPCClustersRow ddd,*clrow=&ddd;
e7de6656 124 clrow = new AliTPCClustersRow("AliTPCcluster");
d7cbafb5 125 cTree->Branch("Segment","AliTPCClustersRow",&clrow,32000,200);
afc42102 126
d7cbafb5 127 const Int_t kMAXZ=fPar->GetMaxTBin()+2;
afc42102 128
73042f01 129 Int_t nclusters=0;
130
c630aafd 131 Int_t nentries = (Int_t)dTree->GetEntries();
132 for (Int_t n=0; n<nentries; n++) {
88cb7938 133
73042f01 134 Int_t sec, row;
c630aafd 135 dTree->GetEvent(n);
88cb7938 136
d7cbafb5 137 if (!fPar->AdjustSectorRow(digarr.GetID(),sec,row)) {
c630aafd 138 Error("Digits2Clusters","!nvalid segment ID ! %d",digarr.GetID());
73042f01 139 continue;
140 }
141
d7cbafb5 142 clrow->SetID(digarr.GetID());
73042f01 143
d7cbafb5 144 cTree->GetBranch("Segment")->SetAddress(&clrow);
145
146 Float_t rx=fPar->GetPadRowRadii(sec,row);
73042f01 147
148 Int_t npads, sign;
149 {
d7cbafb5 150 const Int_t kNIS=fPar->GetNInnerSector(), kNOS=fPar->GetNOuterSector();
73042f01 151 if (sec < kNIS) {
d7cbafb5 152 npads = fPar->GetNPadsLow(row);
73042f01 153 sign = (sec < kNIS/2) ? 1 : -1;
154 } else {
d7cbafb5 155 npads = fPar->GetNPadsUp(row);
73042f01 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();
d7cbafb5 169 if (dig<=fPar->GetZeroSup()) continue;
73042f01 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();
d7cbafb5 217 Float_t w=fPar->GetPadPitchWidth(sec);
73042f01 218 c.SetSigmaY2((s2 + 1./12.)*w*w);
219 if (s2 != 0.) {
220 c.SetSigmaY2(c.GetSigmaY2()*0.108);
d7cbafb5 221 if (sec<fPar->GetNInnerSector()) c.SetSigmaY2(c.GetSigmaY2()*2.07);
73042f01 222 }
223
224 s2 = c.GetSigmaZ2()/c.GetQ() - c.GetZ()*c.GetZ();
d7cbafb5 225 w=fPar->GetZWidth();
73042f01 226 c.SetSigmaZ2((s2 + 1./12.)*w*w);
227 if (s2 != 0.) {
228 c.SetSigmaZ2(c.GetSigmaZ2()*0.169);
d7cbafb5 229 if (sec<fPar->GetNInnerSector()) c.SetSigmaZ2(c.GetSigmaZ2()*1.77);
73042f01 230 }
231
d7cbafb5 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
a1ec4d07 235 c.SetZ(sign*(fPar->GetZLength(sec) - c.GetZ()));
73042f01 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 }
d7cbafb5 253 cTree->Fill();
254
e7de6656 255 clrow->GetArray()->Clear();
73042f01 256
257 nclusters+=ncl;
258
259 delete[] bins;
260 }
261
c630aafd 262 Info("Digits2Cluster","Number of found clusters : %d",nclusters);
f38c8ae5 263
c630aafd 264 return 0;
73042f01 265}
266