Correct sign for calculated b_y (outside measured region).
[u/mrichter/AliRoot.git] / TPC / AliTPCclusterer.cxx
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 /* $Id$ */
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 "AliDigits.h"
29 #include "AliSimDigits.h"
30 #include "AliTPCParam.h"
31 #include "AliTPCClustersRow.h"
32 #include <TTree.h>
33
34 ClassImp(AliTPCclusterer)
35
36
37 //____________________________________________________
38 AliTPCclusterer::AliTPCclusterer(const AliTPCclusterer &param)
39                 :TObject(),
40                  fPar(0)
41 {
42   //
43   // dummy
44   //
45   fPar = param.fPar;
46 }
47 //-----------------------------------------------------
48 AliTPCclusterer & AliTPCclusterer::operator =(const AliTPCclusterer & param)
49 {
50   //
51   // assignment operator - dummy
52   //
53     fPar = param.fPar;
54   return (*this);
55 }
56 //____________________________________________________
57 void AliTPCclusterer::FindPeaks(Int_t k,Int_t max,
58 AliBin *b,Int_t *idx,UInt_t *msk,Int_t& n) {
59   //find local maxima
60   if (n<31)
61   if (IsMaximum(k,max,b)) {
62     idx[n]=k; msk[n]=(2<<n);
63     n++;
64   }
65   b[k].SetMask(0);
66   if (b[k-max].GetMask()&1) FindPeaks(k-max,max,b,idx,msk,n);
67   if (b[k-1  ].GetMask()&1) FindPeaks(k-1  ,max,b,idx,msk,n);
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 }
71
72 void AliTPCclusterer::MarkPeak(Int_t k, Int_t max, AliBin *bins, UInt_t m) {
73   //mark this peak
74   UShort_t q=bins[k].GetQ();
75
76   bins[k].SetMask(bins[k].GetMask()|m); 
77
78   if (bins[k-max].GetQ() <= q)
79      if ((bins[k-max].GetMask()&m) == 0) MarkPeak(k-max,max,bins,m);
80   if (bins[k-1  ].GetQ() <= q)
81      if ((bins[k-1  ].GetMask()&m) == 0) MarkPeak(k-1  ,max,bins,m);
82   if (bins[k+max].GetQ() <= q)
83      if ((bins[k+max].GetMask()&m) == 0) MarkPeak(k+max,max,bins,m);
84   if (bins[k+1  ].GetQ() <= q)
85      if ((bins[k+1  ].GetMask()&m) == 0) MarkPeak(k+1  ,max,bins,m);
86 }
87
88 void AliTPCclusterer::MakeCluster(Int_t k,Int_t max,AliBin *bins,UInt_t m,
89 AliTPCcluster &c) {
90   //make cluster using digits of this peak
91   Float_t q=(Float_t)bins[k].GetQ();
92   Int_t i=k/max, j=k-i*max;
93
94   c.SetQ(c.GetQ()+q);
95   c.SetY(c.GetY()+i*q); 
96   c.SetZ(c.GetZ()+j*q); 
97   c.SetSigmaY2(c.GetSigmaY2()+i*i*q);
98   c.SetSigmaZ2(c.GetSigmaZ2()+j*j*q);
99
100   bins[k].SetMask(0xFFFFFFFE);
101   
102   if (bins[k-max].GetMask() == m) MakeCluster(k-max,max,bins,m,c);
103   if (bins[k-1  ].GetMask() == m) MakeCluster(k-1  ,max,bins,m,c);
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 }
107
108 //_____________________________________________________________________________
109 Int_t AliTPCclusterer::Digits2Clusters(TTree *dTree, TTree *cTree) {
110   //-----------------------------------------------------------------
111   // This is a simple cluster finder.
112   //-----------------------------------------------------------------
113   TBranch *branch=dTree->GetBranch("Segment");
114   if (!branch) {
115      Error("Digits2Cluster","Can't get the branch !");
116      return 1;
117   }
118   AliSimDigits digarr, *dummy=&digarr;
119   branch->SetAddress(&dummy);
120   
121   AliTPCClustersRow ddd,*clrow=&ddd;
122   clrow->SetClass("AliTPCcluster"); clrow->SetArray(1);
123   cTree->Branch("Segment","AliTPCClustersRow",&clrow,32000,200);    
124
125   const Int_t kMAXZ=fPar->GetMaxTBin()+2;
126
127   Int_t nclusters=0;
128
129   Int_t nentries = (Int_t)dTree->GetEntries();
130   for (Int_t n=0; n<nentries; n++) {
131    
132     Int_t sec, row;
133     dTree->GetEvent(n);
134
135     if (!fPar->AdjustSectorRow(digarr.GetID(),sec,row)) {
136        Error("Digits2Clusters","!nvalid segment ID ! %d",digarr.GetID());
137        continue;
138     }
139
140     clrow=new AliTPCClustersRow();
141
142     clrow->SetClass("AliTPCcluster"); clrow->SetArray(1);
143     clrow->SetID(digarr.GetID());
144
145     cTree->GetBranch("Segment")->SetAddress(&clrow);
146
147     Float_t rx=fPar->GetPadRowRadii(sec,row);
148
149     Int_t npads, sign;
150     {
151        const Int_t kNIS=fPar->GetNInnerSector(), kNOS=fPar->GetNOuterSector();
152        if (sec < kNIS) {
153           npads = fPar->GetNPadsLow(row);
154           sign = (sec < kNIS/2) ? 1 : -1;
155        } else {
156           npads = fPar->GetNPadsUp(row);
157           sign = ((sec-kNIS) < kNOS/2) ? 1 : -1;
158        }
159     }
160
161     const Int_t kMAXBIN=kMAXZ*(npads+2);
162     AliBin *bins=new AliBin[kMAXBIN];
163     for (Int_t ii=0;ii<kMAXBIN;ii++) {
164        bins[ii].SetQ(0); bins[ii].SetMask(0xFFFFFFFE);
165     }
166
167     digarr.First();
168     do {
169        Short_t dig=digarr.CurrentDigit();
170        if (dig<=fPar->GetZeroSup()) continue;
171        Int_t j=digarr.CurrentRow()+1, i=digarr.CurrentColumn()+1;
172        bins[i*kMAXZ+j].SetQ(dig);
173        bins[i*kMAXZ+j].SetMask(1);
174     } while (digarr.Next());
175
176     Int_t ncl=0;
177     for (Int_t i=0; i<kMAXBIN; i++) {
178       if ((bins[i].GetMask()&1) == 0) continue;
179       Int_t idx[32]; UInt_t msk[32]; Int_t npeaks=0;
180       FindPeaks(i, kMAXZ, bins, idx, msk, npeaks);
181
182       if (npeaks>30) continue;
183
184       Int_t k,l;
185       for (k=0; k<npeaks-1; k++){//mark adjacent peaks
186         if (idx[k] < 0) continue; //this peak is already removed
187         for (l=k+1; l<npeaks; l++) {
188            if (idx[l] < 0) continue; //this peak is already removed
189            Int_t ki=idx[k]/kMAXZ, kj=idx[k] - ki*kMAXZ;
190            Int_t li=idx[l]/kMAXZ, lj=idx[l] - li*kMAXZ;
191            Int_t di=TMath::Abs(ki - li);
192            Int_t dj=TMath::Abs(kj - lj);
193            if (di>1 || dj>1) continue;
194            if (bins[idx[k]].GetQ() > bins[idx[l]].GetQ()) {
195               msk[l]=msk[k];
196               idx[l]*=-1;
197            } else {
198               msk[k]=msk[l];
199               idx[k]*=-1;
200               break;
201            } 
202         }
203       }
204
205       for (k=0; k<npeaks; k++) {
206         MarkPeak(TMath::Abs(idx[k]), kMAXZ, bins, msk[k]);
207       }
208         
209       for (k=0; k<npeaks; k++) {
210          if (idx[k] < 0) continue; //removed peak
211          AliTPCcluster c;
212          MakeCluster(idx[k], kMAXZ, bins, msk[k], c);
213          if (c.GetQ() < 5) continue; //noise cluster
214          c.SetY(c.GetY()/c.GetQ());
215          c.SetZ(c.GetZ()/c.GetQ());
216
217          Float_t s2 = c.GetSigmaY2()/c.GetQ() - c.GetY()*c.GetY();
218          Float_t w=fPar->GetPadPitchWidth(sec);
219          c.SetSigmaY2((s2 + 1./12.)*w*w);
220          if (s2 != 0.) {
221            c.SetSigmaY2(c.GetSigmaY2()*0.108);
222            if (sec<fPar->GetNInnerSector()) c.SetSigmaY2(c.GetSigmaY2()*2.07);
223          }
224
225          s2 = c.GetSigmaZ2()/c.GetQ() - c.GetZ()*c.GetZ();
226          w=fPar->GetZWidth();
227          c.SetSigmaZ2((s2 + 1./12.)*w*w);
228          if (s2 != 0.) {
229            c.SetSigmaZ2(c.GetSigmaZ2()*0.169);
230            if (sec<fPar->GetNInnerSector()) c.SetSigmaZ2(c.GetSigmaZ2()*1.77);
231          }
232
233          c.SetY((c.GetY() - 0.5 - 0.5*npads)*fPar->GetPadPitchWidth(sec));
234          c.SetZ(fPar->GetZWidth()*(c.GetZ()-1)); 
235          c.SetZ(c.GetZ() - 3.*fPar->GetZSigma()); // PASA delay 
236          c.SetZ(sign*(fPar->GetZLength(sec) - c.GetZ()));
237
238          if (rx<230./250.*TMath::Abs(c.GetZ())) continue;
239
240          Int_t ki=idx[k]/kMAXZ, kj=idx[k] - ki*kMAXZ;
241          c.SetLabel(digarr.GetTrackID(kj-1,ki-1,0),0);
242          c.SetLabel(digarr.GetTrackID(kj-1,ki-1,1),1);
243          c.SetLabel(digarr.GetTrackID(kj-1,ki-1,2),2);
244
245          c.SetQ(bins[idx[k]].GetQ());
246
247          if (ki==1 || ki==npads || kj==1 || kj==kMAXZ-2) {
248            c.SetSigmaY2(c.GetSigmaY2()*25.);
249            c.SetSigmaZ2(c.GetSigmaZ2()*4.);
250          }
251          clrow->InsertCluster(&c); ncl++;
252       }
253     }
254     cTree->Fill();
255
256     delete clrow;
257
258     nclusters+=ncl;
259
260     delete[] bins;  
261   }
262
263   Info("Digits2Cluster","Number of found clusters : %d",nclusters);
264
265   return 0;
266 }
267