1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
18 Revision 1.6 2001/08/30 09:28:48 hristov
19 TTree names are explicitly set via SetName(name) and then Write() is called
21 Revision 1.5 2001/07/20 14:32:44 kowal2
22 Processing of many events possible now
24 Revision 1.4 2001/04/17 08:06:27 hristov
25 Possibility to define the magnetic field in the reconstruction (Yu.Belikov)
27 Revision 1.3 2000/10/05 16:14:01 kowal2
30 Revision 1.2 2000/06/30 12:07:50 kowal2
31 Updated from the TPC-PreRelease branch
33 Revision 1.1.2.1 2000/06/25 08:53:55 kowal2
34 Splitted from AliTPCtracking
38 //-------------------------------------------------------
39 // Implementation of the TPC clusterer
41 // Origin: Jouri Belikov, CERN, Jouri.Belikov@cern.ch
42 //-------------------------------------------------------
44 #include "AliTPCclusterer.h"
45 #include "AliTPCcluster.h"
46 #include <TObjArray.h>
48 #include "AliTPCClustersArray.h"
49 #include "AliTPCClustersRow.h"
50 #include "AliDigits.h"
51 #include "AliSimDigits.h"
52 #include "AliTPCParam.h"
56 void AliTPCclusterer::FindPeaks(Int_t k,Int_t max,
57 AliBin *b,Int_t *idx,UInt_t *msk,Int_t& n) {
60 if (IsMaximum(k,max,b)) {
61 idx[n]=k; msk[n]=(2<<n);
65 if (b[k-max].GetMask()&1) FindPeaks(k-max,max,b,idx,msk,n);
66 if (b[k-1 ].GetMask()&1) FindPeaks(k-1 ,max,b,idx,msk,n);
67 if (b[k+max].GetMask()&1) FindPeaks(k+max,max,b,idx,msk,n);
68 if (b[k+1 ].GetMask()&1) FindPeaks(k+1 ,max,b,idx,msk,n);
71 void AliTPCclusterer::MarkPeak(Int_t k, Int_t max, AliBin *bins, UInt_t m) {
73 UShort_t q=bins[k].GetQ();
75 bins[k].SetMask(bins[k].GetMask()|m);
77 if (bins[k-max].GetQ() <= q)
78 if ((bins[k-max].GetMask()&m) == 0) MarkPeak(k-max,max,bins,m);
79 if (bins[k-1 ].GetQ() <= q)
80 if ((bins[k-1 ].GetMask()&m) == 0) MarkPeak(k-1 ,max,bins,m);
81 if (bins[k+max].GetQ() <= q)
82 if ((bins[k+max].GetMask()&m) == 0) MarkPeak(k+max,max,bins,m);
83 if (bins[k+1 ].GetQ() <= q)
84 if ((bins[k+1 ].GetMask()&m) == 0) MarkPeak(k+1 ,max,bins,m);
87 void AliTPCclusterer::MakeCluster(Int_t k,Int_t max,AliBin *bins,UInt_t m,
89 //make cluster using digits of this peak
90 Float_t q=(Float_t)bins[k].GetQ();
91 Int_t i=k/max, j=k-i*max;
96 c.SetSigmaY2(c.GetSigmaY2()+i*i*q);
97 c.SetSigmaZ2(c.GetSigmaZ2()+j*j*q);
99 bins[k].SetMask(0xFFFFFFFE);
101 if (bins[k-max].GetMask() == m) MakeCluster(k-max,max,bins,m,c);
102 if (bins[k-1 ].GetMask() == m) MakeCluster(k-1 ,max,bins,m,c);
103 if (bins[k+max].GetMask() == m) MakeCluster(k+max,max,bins,m,c);
104 if (bins[k+1 ].GetMask() == m) MakeCluster(k+1 ,max,bins,m,c);
107 //_____________________________________________________________________________
108 void AliTPCclusterer::Digits2Clusters(const AliTPCParam *par, TFile *of, Int_t eventn)
110 //-----------------------------------------------------------------
111 // This is a simple cluster finder.
112 //-----------------------------------------------------------------
113 TDirectory *savedir=gDirectory;
116 cerr<<"AliTPC::Digits2Clusters(): output file not open !\n";
120 const Int_t kMAXZ=par->GetMaxTBin()+2;
126 // for backward compatibility
128 sprintf(dname,"TreeD_75x40_100x60");
129 sprintf(cname,"TreeC_TPC");
132 sprintf(dname,"TreeD_75x40_100x60_%d",eventn);
133 sprintf(cname,"TreeC_TPC_%d",eventn);
135 TTree *t = (TTree *)gDirectory->Get(dname);
137 AliSimDigits digarr, *dummy=&digarr;
138 t->GetBranch("Segment")->SetAddress(&dummy);
139 Stat_t nentries = t->GetEntries();
143 ((AliTPCParam*)par)->Write(par->GetTitle());
144 AliTPCClustersArray carray;
146 carray.SetClusterType("AliTPCcluster");
153 for (Int_t n=0; n<nentries; n++) {
156 if (!par->AdjustSectorRow(digarr.GetID(),sec,row)) {
157 cerr<<"AliTPC warning: invalid segment ID ! "<<digarr.GetID()<<endl;
161 AliTPCClustersRow *clrow=carray.CreateRow(sec,row);
163 Float_t rx=par->GetPadRowRadii(sec,row);
167 const Int_t kNIS=par->GetNInnerSector(), kNOS=par->GetNOuterSector();
169 npads = par->GetNPadsLow(row);
170 sign = (sec < kNIS/2) ? 1 : -1;
172 npads = par->GetNPadsUp(row);
173 sign = ((sec-kNIS) < kNOS/2) ? 1 : -1;
177 const Int_t kMAXBIN=kMAXZ*(npads+2);
178 AliBin *bins=new AliBin[kMAXBIN];
179 for (Int_t ii=0;ii<kMAXBIN;ii++) {
180 bins[ii].SetQ(0); bins[ii].SetMask(0xFFFFFFFE);
185 Short_t dig=digarr.CurrentDigit();
186 if (dig<=par->GetZeroSup()) continue;
187 Int_t j=digarr.CurrentRow()+1, i=digarr.CurrentColumn()+1;
188 bins[i*kMAXZ+j].SetQ(dig);
189 bins[i*kMAXZ+j].SetMask(1);
190 } while (digarr.Next());
193 for (Int_t i=0; i<kMAXBIN; i++) {
194 if ((bins[i].GetMask()&1) == 0) continue;
195 Int_t idx[32]; UInt_t msk[32]; Int_t npeaks=0;
196 FindPeaks(i, kMAXZ, bins, idx, msk, npeaks);
198 if (npeaks>30) continue;
201 for (k=0; k<npeaks-1; k++){//mark adjacent peaks
202 if (idx[k] < 0) continue; //this peak is already removed
203 for (l=k+1; l<npeaks; l++) {
204 if (idx[l] < 0) continue; //this peak is already removed
205 Int_t ki=idx[k]/kMAXZ, kj=idx[k] - ki*kMAXZ;
206 Int_t li=idx[l]/kMAXZ, lj=idx[l] - li*kMAXZ;
207 Int_t di=TMath::Abs(ki - li);
208 Int_t dj=TMath::Abs(kj - lj);
209 if (di>1 || dj>1) continue;
210 if (bins[idx[k]].GetQ() > bins[idx[l]].GetQ()) {
221 for (k=0; k<npeaks; k++) {
222 MarkPeak(TMath::Abs(idx[k]), kMAXZ, bins, msk[k]);
225 for (k=0; k<npeaks; k++) {
226 if (idx[k] < 0) continue; //removed peak
228 MakeCluster(idx[k], kMAXZ, bins, msk[k], c);
229 if (c.GetQ() < 5) continue; //noise cluster
230 c.SetY(c.GetY()/c.GetQ());
231 c.SetZ(c.GetZ()/c.GetQ());
233 Float_t s2 = c.GetSigmaY2()/c.GetQ() - c.GetY()*c.GetY();
234 Float_t w=par->GetPadPitchWidth(sec);
235 c.SetSigmaY2((s2 + 1./12.)*w*w);
237 c.SetSigmaY2(c.GetSigmaY2()*0.108);
238 if (sec<par->GetNInnerSector()) c.SetSigmaY2(c.GetSigmaY2()*2.07);
241 s2 = c.GetSigmaZ2()/c.GetQ() - c.GetZ()*c.GetZ();
243 c.SetSigmaZ2((s2 + 1./12.)*w*w);
245 c.SetSigmaZ2(c.GetSigmaZ2()*0.169);
246 if (sec<par->GetNInnerSector()) c.SetSigmaZ2(c.GetSigmaZ2()*1.77);
249 c.SetY((c.GetY() - 0.5 - 0.5*npads)*par->GetPadPitchWidth(sec));
250 c.SetZ(par->GetZWidth()*(c.GetZ()-1));
251 c.SetZ(c.GetZ() - 3.*par->GetZSigma()); // PASA delay
252 c.SetZ(sign*(par->GetZLength() - c.GetZ()));
254 if (rx<230./250.*TMath::Abs(c.GetZ())) continue;
256 Int_t ki=idx[k]/kMAXZ, kj=idx[k] - ki*kMAXZ;
257 c.SetLabel(digarr.GetTrackID(kj-1,ki-1,0),0);
258 c.SetLabel(digarr.GetTrackID(kj-1,ki-1,1),1);
259 c.SetLabel(digarr.GetTrackID(kj-1,ki-1,2),2);
261 c.SetQ(bins[idx[k]].GetQ());
263 if (ki==1 || ki==npads || kj==1 || kj==kMAXZ-2) {
264 c.SetSigmaY2(c.GetSigmaY2()*25.);
265 c.SetSigmaZ2(c.GetSigmaZ2()*4.);
267 clrow->InsertCluster(&c); ncl++;
270 carray.StoreRow(sec,row);
271 carray.ClearRow(sec,row);
273 //cerr<<"sector, row, compressed digits, clusters: "
274 //<<sec<<' '<<row<<' '<<digarr.GetSize()<<' '<<ncl<<" \r";
281 cerr<<"Number of found clusters : "<<nclusters<<" \n";
283 carray.GetTree()->SetName(cname);
284 carray.GetTree()->Write();
287 delete t; //Thanks to Mariana Bondila