Protection against nonexisting input tree
[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
16/*
17$Log$
6e4d905e 18Revision 1.8 2002/03/29 06:57:45 kowal2
19Restored backward compatibility to use the hits from Dec. 2000 production.
20
7a09f434 21Revision 1.7 2001/10/21 19:04:55 hristov
22Several patches were done to adapt the barel reconstruction to the multi-event case. Some memory leaks were corrected. (Yu.Belikov)
23
f38c8ae5 24Revision 1.6 2001/08/30 09:28:48 hristov
25TTree names are explicitly set via SetName(name) and then Write() is called
26
6bc3458d 27Revision 1.5 2001/07/20 14:32:44 kowal2
28Processing of many events possible now
29
afc42102 30Revision 1.4 2001/04/17 08:06:27 hristov
31Possibility to define the magnetic field in the reconstruction (Yu.Belikov)
32
9b280d80 33Revision 1.3 2000/10/05 16:14:01 kowal2
34Forward declarations.
35
c3c63118 36Revision 1.2 2000/06/30 12:07:50 kowal2
37Updated from the TPC-PreRelease branch
38
73042f01 39Revision 1.1.2.1 2000/06/25 08:53:55 kowal2
40Splitted from AliTPCtracking
41
42*/
43
44//-------------------------------------------------------
45// Implementation of the TPC clusterer
46//
47// Origin: Jouri Belikov, CERN, Jouri.Belikov@cern.ch
48//-------------------------------------------------------
49
50#include "AliTPCclusterer.h"
51#include "AliTPCcluster.h"
52#include <TObjArray.h>
53#include <TFile.h>
54#include "AliTPCClustersArray.h"
55#include "AliTPCClustersRow.h"
56#include "AliDigits.h"
57#include "AliSimDigits.h"
58#include "AliTPCParam.h"
59#include <iostream.h>
c3c63118 60#include <TTree.h>
73042f01 61
62void AliTPCclusterer::FindPeaks(Int_t k,Int_t max,
63AliBin *b,Int_t *idx,UInt_t *msk,Int_t& n) {
64 //find local maxima
65 if (n<31)
66 if (IsMaximum(k,max,b)) {
67 idx[n]=k; msk[n]=(2<<n);
68 n++;
69 }
70 b[k].SetMask(0);
71 if (b[k-max].GetMask()&1) FindPeaks(k-max,max,b,idx,msk,n);
72 if (b[k-1 ].GetMask()&1) FindPeaks(k-1 ,max,b,idx,msk,n);
73 if (b[k+max].GetMask()&1) FindPeaks(k+max,max,b,idx,msk,n);
74 if (b[k+1 ].GetMask()&1) FindPeaks(k+1 ,max,b,idx,msk,n);
75}
76
77void AliTPCclusterer::MarkPeak(Int_t k, Int_t max, AliBin *bins, UInt_t m) {
78 //mark this peak
79 UShort_t q=bins[k].GetQ();
80
81 bins[k].SetMask(bins[k].GetMask()|m);
82
83 if (bins[k-max].GetQ() <= q)
84 if ((bins[k-max].GetMask()&m) == 0) MarkPeak(k-max,max,bins,m);
85 if (bins[k-1 ].GetQ() <= q)
86 if ((bins[k-1 ].GetMask()&m) == 0) MarkPeak(k-1 ,max,bins,m);
87 if (bins[k+max].GetQ() <= q)
88 if ((bins[k+max].GetMask()&m) == 0) MarkPeak(k+max,max,bins,m);
89 if (bins[k+1 ].GetQ() <= q)
90 if ((bins[k+1 ].GetMask()&m) == 0) MarkPeak(k+1 ,max,bins,m);
91}
92
93void AliTPCclusterer::MakeCluster(Int_t k,Int_t max,AliBin *bins,UInt_t m,
94AliTPCcluster &c) {
95 //make cluster using digits of this peak
96 Float_t q=(Float_t)bins[k].GetQ();
97 Int_t i=k/max, j=k-i*max;
98
99 c.SetQ(c.GetQ()+q);
100 c.SetY(c.GetY()+i*q);
101 c.SetZ(c.GetZ()+j*q);
102 c.SetSigmaY2(c.GetSigmaY2()+i*i*q);
103 c.SetSigmaZ2(c.GetSigmaZ2()+j*j*q);
104
105 bins[k].SetMask(0xFFFFFFFE);
106
107 if (bins[k-max].GetMask() == m) MakeCluster(k-max,max,bins,m,c);
108 if (bins[k-1 ].GetMask() == m) MakeCluster(k-1 ,max,bins,m,c);
109 if (bins[k+max].GetMask() == m) MakeCluster(k+max,max,bins,m,c);
110 if (bins[k+1 ].GetMask() == m) MakeCluster(k+1 ,max,bins,m,c);
111}
112
113//_____________________________________________________________________________
afc42102 114void AliTPCclusterer::Digits2Clusters(const AliTPCParam *par, TFile *of, Int_t eventn)
73042f01 115{
116 //-----------------------------------------------------------------
117 // This is a simple cluster finder.
118 //-----------------------------------------------------------------
119 TDirectory *savedir=gDirectory;
120
121 if (!of->IsOpen()) {
122 cerr<<"AliTPC::Digits2Clusters(): output file not open !\n";
123 return;
124 }
125
126 const Int_t kMAXZ=par->GetMaxTBin()+2;
127
afc42102 128 char dname[100];
129 char cname[100];
130 if (eventn==-1) {
131
132 // for backward compatibility
133
7a09f434 134 sprintf(dname,"TreeD_75x40_100x60_150x60");
afc42102 135 sprintf(cname,"TreeC_TPC");
136 }
137 else {
7a09f434 138 sprintf(dname,"TreeD_75x40_100x60_150x60_%d",eventn);
afc42102 139 sprintf(cname,"TreeC_TPC_%d",eventn);
140 }
141 TTree *t = (TTree *)gDirectory->Get(dname);
142
6e4d905e 143 if (!t) {
144 cerr<<"Input tree with "<<dname<<" not found"<<endl;
145 return;
146 }
147
73042f01 148 AliSimDigits digarr, *dummy=&digarr;
149 t->GetBranch("Segment")->SetAddress(&dummy);
150 Stat_t nentries = t->GetEntries();
151
152 of->cd();
153
154 ((AliTPCParam*)par)->Write(par->GetTitle());
155 AliTPCClustersArray carray;
156 carray.Setup(par);
157 carray.SetClusterType("AliTPCcluster");
158 carray.MakeTree();
159
afc42102 160
161
73042f01 162 Int_t nclusters=0;
163
164 for (Int_t n=0; n<nentries; n++) {
165 t->GetEvent(n);
166 Int_t sec, row;
167 if (!par->AdjustSectorRow(digarr.GetID(),sec,row)) {
168 cerr<<"AliTPC warning: invalid segment ID ! "<<digarr.GetID()<<endl;
169 continue;
170 }
171
172 AliTPCClustersRow *clrow=carray.CreateRow(sec,row);
173
174 Float_t rx=par->GetPadRowRadii(sec,row);
175
176 Int_t npads, sign;
177 {
178 const Int_t kNIS=par->GetNInnerSector(), kNOS=par->GetNOuterSector();
179 if (sec < kNIS) {
180 npads = par->GetNPadsLow(row);
181 sign = (sec < kNIS/2) ? 1 : -1;
182 } else {
183 npads = par->GetNPadsUp(row);
184 sign = ((sec-kNIS) < kNOS/2) ? 1 : -1;
185 }
186 }
187
188 const Int_t kMAXBIN=kMAXZ*(npads+2);
189 AliBin *bins=new AliBin[kMAXBIN];
190 for (Int_t ii=0;ii<kMAXBIN;ii++) {
191 bins[ii].SetQ(0); bins[ii].SetMask(0xFFFFFFFE);
192 }
193
194 digarr.First();
195 do {
196 Short_t dig=digarr.CurrentDigit();
197 if (dig<=par->GetZeroSup()) continue;
198 Int_t j=digarr.CurrentRow()+1, i=digarr.CurrentColumn()+1;
199 bins[i*kMAXZ+j].SetQ(dig);
200 bins[i*kMAXZ+j].SetMask(1);
201 } while (digarr.Next());
202
203 Int_t ncl=0;
204 for (Int_t i=0; i<kMAXBIN; i++) {
205 if ((bins[i].GetMask()&1) == 0) continue;
206 Int_t idx[32]; UInt_t msk[32]; Int_t npeaks=0;
207 FindPeaks(i, kMAXZ, bins, idx, msk, npeaks);
208
209 if (npeaks>30) continue;
210
211 Int_t k,l;
212 for (k=0; k<npeaks-1; k++){//mark adjacent peaks
213 if (idx[k] < 0) continue; //this peak is already removed
214 for (l=k+1; l<npeaks; l++) {
215 if (idx[l] < 0) continue; //this peak is already removed
216 Int_t ki=idx[k]/kMAXZ, kj=idx[k] - ki*kMAXZ;
217 Int_t li=idx[l]/kMAXZ, lj=idx[l] - li*kMAXZ;
218 Int_t di=TMath::Abs(ki - li);
219 Int_t dj=TMath::Abs(kj - lj);
220 if (di>1 || dj>1) continue;
221 if (bins[idx[k]].GetQ() > bins[idx[l]].GetQ()) {
222 msk[l]=msk[k];
223 idx[l]*=-1;
224 } else {
225 msk[k]=msk[l];
226 idx[k]*=-1;
227 break;
228 }
229 }
230 }
231
232 for (k=0; k<npeaks; k++) {
233 MarkPeak(TMath::Abs(idx[k]), kMAXZ, bins, msk[k]);
234 }
235
236 for (k=0; k<npeaks; k++) {
237 if (idx[k] < 0) continue; //removed peak
238 AliTPCcluster c;
239 MakeCluster(idx[k], kMAXZ, bins, msk[k], c);
240 if (c.GetQ() < 5) continue; //noise cluster
241 c.SetY(c.GetY()/c.GetQ());
242 c.SetZ(c.GetZ()/c.GetQ());
243
244 Float_t s2 = c.GetSigmaY2()/c.GetQ() - c.GetY()*c.GetY();
245 Float_t w=par->GetPadPitchWidth(sec);
246 c.SetSigmaY2((s2 + 1./12.)*w*w);
247 if (s2 != 0.) {
248 c.SetSigmaY2(c.GetSigmaY2()*0.108);
249 if (sec<par->GetNInnerSector()) c.SetSigmaY2(c.GetSigmaY2()*2.07);
250 }
251
252 s2 = c.GetSigmaZ2()/c.GetQ() - c.GetZ()*c.GetZ();
253 w=par->GetZWidth();
254 c.SetSigmaZ2((s2 + 1./12.)*w*w);
255 if (s2 != 0.) {
256 c.SetSigmaZ2(c.GetSigmaZ2()*0.169);
257 if (sec<par->GetNInnerSector()) c.SetSigmaZ2(c.GetSigmaZ2()*1.77);
258 }
259
260 c.SetY((c.GetY() - 0.5 - 0.5*npads)*par->GetPadPitchWidth(sec));
261 c.SetZ(par->GetZWidth()*(c.GetZ()-1));
262 c.SetZ(c.GetZ() - 3.*par->GetZSigma()); // PASA delay
263 c.SetZ(sign*(par->GetZLength() - c.GetZ()));
264
265 if (rx<230./250.*TMath::Abs(c.GetZ())) continue;
266
267 Int_t ki=idx[k]/kMAXZ, kj=idx[k] - ki*kMAXZ;
268 c.SetLabel(digarr.GetTrackID(kj-1,ki-1,0),0);
269 c.SetLabel(digarr.GetTrackID(kj-1,ki-1,1),1);
270 c.SetLabel(digarr.GetTrackID(kj-1,ki-1,2),2);
271
272 c.SetQ(bins[idx[k]].GetQ());
273
274 if (ki==1 || ki==npads || kj==1 || kj==kMAXZ-2) {
275 c.SetSigmaY2(c.GetSigmaY2()*25.);
276 c.SetSigmaZ2(c.GetSigmaZ2()*4.);
277 }
278 clrow->InsertCluster(&c); ncl++;
279 }
280 }
281 carray.StoreRow(sec,row);
282 carray.ClearRow(sec,row);
283
9b280d80 284 //cerr<<"sector, row, compressed digits, clusters: "
285 //<<sec<<' '<<row<<' '<<digarr.GetSize()<<' '<<ncl<<" \r";
73042f01 286
287 nclusters+=ncl;
288
289 delete[] bins;
290 }
291
292 cerr<<"Number of found clusters : "<<nclusters<<" \n";
293
6bc3458d 294 carray.GetTree()->SetName(cname);
295 carray.GetTree()->Write();
73042f01 296 savedir->cd();
f38c8ae5 297
298 delete t; //Thanks to Mariana Bondila
73042f01 299}
300