Transition to NewIO
[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>
28#include "AliTPCClustersArray.h"
29#include "AliTPCClustersRow.h"
30#include "AliDigits.h"
31#include "AliSimDigits.h"
32#include "AliTPCParam.h"
19364939 33#include <Riostream.h>
c3c63118 34#include <TTree.h>
88cb7938 35#include "AliRunLoader.h"
36#include "AliLoader.h"
73042f01 37
38void AliTPCclusterer::FindPeaks(Int_t k,Int_t max,
39AliBin *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
53void 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
69void AliTPCclusterer::MakeCluster(Int_t k,Int_t max,AliBin *bins,UInt_t m,
70AliTPCcluster &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//_____________________________________________________________________________
88cb7938 90void AliTPCclusterer::Digits2Clusters(const AliTPCParam *par, AliLoader *of, Int_t eventn)
73042f01 91{
92 //-----------------------------------------------------------------
93 // This is a simple cluster finder.
94 //-----------------------------------------------------------------
88cb7938 95 AliRunLoader* rl = (AliRunLoader*)of->GetEventFolder()->FindObject(AliRunLoader::fgkRunLoaderName);
96 rl->GetEvent(eventn);
73042f01 97 TDirectory *savedir=gDirectory;
98
88cb7938 99 if (of->TreeR() == 0x0) of->MakeTree("R");
100
101
102 if (of == 0x0)
103 {
73042f01 104 cerr<<"AliTPC::Digits2Clusters(): output file not open !\n";
105 return;
88cb7938 106 }
73042f01 107
108 const Int_t kMAXZ=par->GetMaxTBin()+2;
109
88cb7938 110
111 TTree *t = (TTree *)of->TreeD();
afc42102 112
6e4d905e 113 if (!t) {
88cb7938 114 cerr<<"Input tree with digits not found"<<endl;
6e4d905e 115 return;
116 }
117
73042f01 118 AliSimDigits digarr, *dummy=&digarr;
119 t->GetBranch("Segment")->SetAddress(&dummy);
120 Stat_t nentries = t->GetEntries();
121
88cb7938 122 cout<<"Got "<<nentries<<" from TreeD"<<endl;
73042f01 123
88cb7938 124// ((AliTPCParam*)par)->Write(par->GetTitle());
125
73042f01 126 AliTPCClustersArray carray;
127 carray.Setup(par);
128 carray.SetClusterType("AliTPCcluster");
88cb7938 129
130 TTree* treeR = of->TreeR();
131 carray.MakeTree(treeR);
73042f01 132
afc42102 133
134
73042f01 135 Int_t nclusters=0;
136
88cb7938 137 for (Int_t n=0; n<nentries; n++)
138 {
139
73042f01 140 Int_t sec, row;
88cb7938 141 t->GetEvent(n);
142
73042f01 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
9b280d80 260 //cerr<<"sector, row, compressed digits, clusters: "
261 //<<sec<<' '<<row<<' '<<digarr.GetSize()<<' '<<ncl<<" \r";
73042f01 262
263 nclusters+=ncl;
264
265 delete[] bins;
266 }
267
268 cerr<<"Number of found clusters : "<<nclusters<<" \n";
269
88cb7938 270 of->WriteRecPoints("OVERWRITE");
271
73042f01 272 savedir->cd();
f38c8ae5 273
88cb7938 274// delete t; //Thanks to Mariana Bondila
73042f01 275}
276