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 //-------------------------------------------------------
19 // Implementation of the TPC clusterer
21 // Origin: Jouri Belikov, CERN, Jouri.Belikov@cern.ch
22 //-------------------------------------------------------
24 #include "AliTPCclusterer.h"
25 #include "AliTPCcluster.h"
26 #include <TObjArray.h>
28 #include "AliDigits.h"
29 #include "AliSimDigits.h"
30 #include "AliTPCParam.h"
31 #include "AliTPCClustersRow.h"
34 ClassImp(AliTPCclusterer)
36 AliTPCclusterer::AliTPCclusterer(const AliTPCParam *par) {
37 //-------------------------------------------------------
38 // The main constructor
39 //-------------------------------------------------------
43 void AliTPCclusterer::FindPeaks(Int_t k,Int_t max,
44 AliBin *b,Int_t *idx,UInt_t *msk,Int_t& n) {
47 if (IsMaximum(k,max,b)) {
48 idx[n]=k; msk[n]=(2<<n);
52 if (b[k-max].GetMask()&1) FindPeaks(k-max,max,b,idx,msk,n);
53 if (b[k-1 ].GetMask()&1) FindPeaks(k-1 ,max,b,idx,msk,n);
54 if (b[k+max].GetMask()&1) FindPeaks(k+max,max,b,idx,msk,n);
55 if (b[k+1 ].GetMask()&1) FindPeaks(k+1 ,max,b,idx,msk,n);
58 void AliTPCclusterer::MarkPeak(Int_t k, Int_t max, AliBin *bins, UInt_t m) {
60 UShort_t q=bins[k].GetQ();
62 bins[k].SetMask(bins[k].GetMask()|m);
64 if (bins[k-max].GetQ() <= q)
65 if ((bins[k-max].GetMask()&m) == 0) MarkPeak(k-max,max,bins,m);
66 if (bins[k-1 ].GetQ() <= q)
67 if ((bins[k-1 ].GetMask()&m) == 0) MarkPeak(k-1 ,max,bins,m);
68 if (bins[k+max].GetQ() <= q)
69 if ((bins[k+max].GetMask()&m) == 0) MarkPeak(k+max,max,bins,m);
70 if (bins[k+1 ].GetQ() <= q)
71 if ((bins[k+1 ].GetMask()&m) == 0) MarkPeak(k+1 ,max,bins,m);
74 void AliTPCclusterer::MakeCluster(Int_t k,Int_t max,AliBin *bins,UInt_t m,
76 //make cluster using digits of this peak
77 Float_t q=(Float_t)bins[k].GetQ();
78 Int_t i=k/max, j=k-i*max;
83 c.SetSigmaY2(c.GetSigmaY2()+i*i*q);
84 c.SetSigmaZ2(c.GetSigmaZ2()+j*j*q);
86 bins[k].SetMask(0xFFFFFFFE);
88 if (bins[k-max].GetMask() == m) MakeCluster(k-max,max,bins,m,c);
89 if (bins[k-1 ].GetMask() == m) MakeCluster(k-1 ,max,bins,m,c);
90 if (bins[k+max].GetMask() == m) MakeCluster(k+max,max,bins,m,c);
91 if (bins[k+1 ].GetMask() == m) MakeCluster(k+1 ,max,bins,m,c);
94 //_____________________________________________________________________________
95 Int_t AliTPCclusterer::Digits2Clusters(TTree *dTree, TTree *cTree) {
96 //-----------------------------------------------------------------
97 // This is a simple cluster finder.
98 //-----------------------------------------------------------------
99 TBranch *branch=dTree->GetBranch("Segment");
101 Error("Digits2Cluster","Can't get the branch !");
104 AliSimDigits digarr, *dummy=&digarr;
105 branch->SetAddress(&dummy);
107 AliTPCClustersRow ddd,*clrow=&ddd;
108 clrow->SetClass("AliTPCcluster"); clrow->SetArray(1);
109 cTree->Branch("Segment","AliTPCClustersRow",&clrow,32000,200);
111 const Int_t kMAXZ=fPar->GetMaxTBin()+2;
115 Int_t nentries = (Int_t)dTree->GetEntries();
116 for (Int_t n=0; n<nentries; n++) {
121 if (!fPar->AdjustSectorRow(digarr.GetID(),sec,row)) {
122 Error("Digits2Clusters","!nvalid segment ID ! %d",digarr.GetID());
126 clrow=new AliTPCClustersRow();
128 clrow->SetClass("AliTPCcluster"); clrow->SetArray(1);
129 clrow->SetID(digarr.GetID());
131 cTree->GetBranch("Segment")->SetAddress(&clrow);
133 Float_t rx=fPar->GetPadRowRadii(sec,row);
137 const Int_t kNIS=fPar->GetNInnerSector(), kNOS=fPar->GetNOuterSector();
139 npads = fPar->GetNPadsLow(row);
140 sign = (sec < kNIS/2) ? 1 : -1;
142 npads = fPar->GetNPadsUp(row);
143 sign = ((sec-kNIS) < kNOS/2) ? 1 : -1;
147 const Int_t kMAXBIN=kMAXZ*(npads+2);
148 AliBin *bins=new AliBin[kMAXBIN];
149 for (Int_t ii=0;ii<kMAXBIN;ii++) {
150 bins[ii].SetQ(0); bins[ii].SetMask(0xFFFFFFFE);
155 Short_t dig=digarr.CurrentDigit();
156 if (dig<=fPar->GetZeroSup()) continue;
157 Int_t j=digarr.CurrentRow()+1, i=digarr.CurrentColumn()+1;
158 bins[i*kMAXZ+j].SetQ(dig);
159 bins[i*kMAXZ+j].SetMask(1);
160 } while (digarr.Next());
163 for (Int_t i=0; i<kMAXBIN; i++) {
164 if ((bins[i].GetMask()&1) == 0) continue;
165 Int_t idx[32]; UInt_t msk[32]; Int_t npeaks=0;
166 FindPeaks(i, kMAXZ, bins, idx, msk, npeaks);
168 if (npeaks>30) continue;
171 for (k=0; k<npeaks-1; k++){//mark adjacent peaks
172 if (idx[k] < 0) continue; //this peak is already removed
173 for (l=k+1; l<npeaks; l++) {
174 if (idx[l] < 0) continue; //this peak is already removed
175 Int_t ki=idx[k]/kMAXZ, kj=idx[k] - ki*kMAXZ;
176 Int_t li=idx[l]/kMAXZ, lj=idx[l] - li*kMAXZ;
177 Int_t di=TMath::Abs(ki - li);
178 Int_t dj=TMath::Abs(kj - lj);
179 if (di>1 || dj>1) continue;
180 if (bins[idx[k]].GetQ() > bins[idx[l]].GetQ()) {
191 for (k=0; k<npeaks; k++) {
192 MarkPeak(TMath::Abs(idx[k]), kMAXZ, bins, msk[k]);
195 for (k=0; k<npeaks; k++) {
196 if (idx[k] < 0) continue; //removed peak
198 MakeCluster(idx[k], kMAXZ, bins, msk[k], c);
199 if (c.GetQ() < 5) continue; //noise cluster
200 c.SetY(c.GetY()/c.GetQ());
201 c.SetZ(c.GetZ()/c.GetQ());
203 Float_t s2 = c.GetSigmaY2()/c.GetQ() - c.GetY()*c.GetY();
204 Float_t w=fPar->GetPadPitchWidth(sec);
205 c.SetSigmaY2((s2 + 1./12.)*w*w);
207 c.SetSigmaY2(c.GetSigmaY2()*0.108);
208 if (sec<fPar->GetNInnerSector()) c.SetSigmaY2(c.GetSigmaY2()*2.07);
211 s2 = c.GetSigmaZ2()/c.GetQ() - c.GetZ()*c.GetZ();
213 c.SetSigmaZ2((s2 + 1./12.)*w*w);
215 c.SetSigmaZ2(c.GetSigmaZ2()*0.169);
216 if (sec<fPar->GetNInnerSector()) c.SetSigmaZ2(c.GetSigmaZ2()*1.77);
219 c.SetY((c.GetY() - 0.5 - 0.5*npads)*fPar->GetPadPitchWidth(sec));
220 c.SetZ(fPar->GetZWidth()*(c.GetZ()-1));
221 c.SetZ(c.GetZ() - 3.*fPar->GetZSigma()); // PASA delay
222 c.SetZ(sign*(fPar->GetZLength() - c.GetZ()));
224 if (rx<230./250.*TMath::Abs(c.GetZ())) continue;
226 Int_t ki=idx[k]/kMAXZ, kj=idx[k] - ki*kMAXZ;
227 c.SetLabel(digarr.GetTrackID(kj-1,ki-1,0),0);
228 c.SetLabel(digarr.GetTrackID(kj-1,ki-1,1),1);
229 c.SetLabel(digarr.GetTrackID(kj-1,ki-1,2),2);
231 c.SetQ(bins[idx[k]].GetQ());
233 if (ki==1 || ki==npads || kj==1 || kj==kMAXZ-2) {
234 c.SetSigmaY2(c.GetSigmaY2()*25.);
235 c.SetSigmaZ2(c.GetSigmaZ2()*4.);
237 clrow->InsertCluster(&c); ncl++;
249 Info("Digits2Cluster","Number of found clusters : %d",nclusters);