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