Comment changes for htlm docu.
[u/mrichter/AliRoot.git] / HLT / misc / AliL3VHDLClusterFinder.cxx
CommitLineData
1e843b89 1// $Id$
2
8e1f5064 3// Author: Constantin Loizides <mailto:loizides@ikf.physik.uni-frankfurt.de>
1e843b89 4//*-- Copyright & copy CL
8e1f5064 5
6#include <math.h>
7#include <time.h>
8#include <iostream.h>
9#include <stdio.h>
10#include <string.h>
11
12#include "AliL3VHDLClusterFinder.h"
13
14/** \class AliL3VHDLClusterFinder
15//<pre>
16//____________________________________________________
17// AliL3VHDLClusterFinder
18//
19// The current VHDL cluster finder for HLT
20// Based on STAR L3
21//
22// Most important parameters:
23// fThreshold - threshold for noise clusters
24// fMatch - length in time for overlapping sequences
6aec85a4 25//</pre>
26*/
8e1f5064 27
28ClassImp(AliL3VHDLClusterFinder)
29
30AliL3VHDLClusterFinder::AliL3VHDLClusterFinder()
31{
32 fMatch = 4;
33 fThreshold = 10;
34 fMinMerge = 1;
35 fNClusters=0;
36
37 fXYErr = 0.2;
38 fZErr = 0.3;
39 fDeconvPad = kTRUE;
40 fDeconvTime = kTRUE;
41 fstdout = kFALSE;
42 fcalcerr = kTRUE;
43
44 Clear();
45#ifdef DEBUG
46 fdeb=fopen("vhdlclusterfinder.debug","w");
47 //fdeb=stderr;
48#endif
49}
50
51AliL3VHDLClusterFinder::~AliL3VHDLClusterFinder()
52{
53#ifdef DEBUG
54 fclose(fdeb);
55#endif
56}
57
58void AliL3VHDLClusterFinder::ProcessDigits()
59{
60 //Loop over data like the analyzer of the VHDL code
61 const UChar_t n=255;
62 UShort_t rrow=0,rtime=0;
63 UChar_t rpad=0,i=n;
64 UShort_t *charges=new UShort_t[n];
65 Int_t tc=0,mp=0,mt=0,sp=0,st=0;
66
67 fNClusters=0;
68 fRow=0;
69 fPad=0;
70 Clear();
71
72 //loop over input data
73 while(fAltromem.ReadSequence(rrow,rpad,rtime,i,&charges)){
74 tc=0;mp=0;mt=0;sp=0;st=0;
75#if 0
76 cout << "Padrow " << (int)rrow << " pad " << (int)rpad << " time " <<(int)rtime << " charges ";
77 for(UChar_t ii=0;ii<i;ii++) cout << (int)charges[ii] << " ";
78 cout << endl;
79#endif
80#ifdef DEBUG
81 fprintf(fdeb,"ProcessDigits: Input Data: %d %d %d charges:",(int)rrow,(int)rpad,(int)rtime);
82 for(UChar_t ii=0;ii<i;ii++) fprintf(fdeb," %d",(int)charges[ii]);
83 fprintf(fdeb,"\n");
84#endif
85
86 fNRow=rrow;
87 fNPad=rpad;
88
89 //calculate sequence values
90 //no deconvulution so far
91 for(UChar_t ii=0;ii<i;ii++){
92 tc+=charges[ii];
93 mt+=(rtime-ii)*charges[ii];
94 st+=(rtime-ii)*(rtime-ii)*charges[ii];
95 }
96 mp=rpad*tc;
97 sp=rpad*rpad*tc;
98
99 fSeq.fTotalCharge=tc;
100 fSeq.fPad=mp;
101 fSeq.fTime=mt;
102 fSeq.fPad2=sp;
103 fSeq.fTime2=st;
104 fSeq.fMean=0;
6aec85a4 105 //if(tc!=0) fSeq.fMean=mt/tc;
106 if(tc!=0) fSeq.fMean=rtime-i/2;
8e1f5064 107 fSeq.fMerge=0;
108 fSeq.fRow=rrow;
109 fSeq.fLastPad=rpad;
110
111 //work on this sequence
112 ProcessSequence();
113 //output one cluster
114 OutputMemory();
115
116#ifdef DEBUG
117 fflush(fdeb);
118#endif
119 i=n; //store size of charges array
120 } //loop over data
121
122 //flush everything left
123 FlushMemory();
124 while(fOP!=fFP) OutputMemory();
125}
126
127void AliL3VHDLClusterFinder::ProcessSequence()
128{
129 if(fNRow!=fRow) FlushMemory();
130 else if(fNPad==fPad+1) PrepareMemory();
131 else if(fNPad!=fPad) FlushMemory();
132
133 //store new row and pad values
134 fRow=fNRow;
135 fPad=fNPad;
136
137#ifdef DEBUG
138 fprintf(fdeb,"ProcessSequence: Mean=%d TC=%d ",fSeq.fMean,fSeq.fTotalCharge);
139#endif
140
141 CompareSeq(); //merge or insert
142}
143
144void AliL3VHDLClusterFinder::PrepareMemory()
145{
146#ifdef DEBUG
147 fprintf(fdeb,"PrepareMemory %d %d %d\n",fRP,fEP,fWP);
148#endif
149 fRP=fEP;
150 fEP=fWP;
151}
152
153void AliL3VHDLClusterFinder::FlushMemory()
154{
155#ifdef DEBUG
156 fprintf(fdeb,"FlushMemory %d %d %d %d\n",fFP,fRP,fEP,fWP);
157#endif
158 fFP=fWP;
159 fRP=fWP;
160 fEP=fWP;
161}
162
163void AliL3VHDLClusterFinder::CompareSeq()
164{
165
166 while(fRP!=fEP){
167 Int_t diff=fSeqs[fPList[fRP]].fMean-fSeq.fMean;
168
169 if(diff>fMatch){ //no match
170 IncRPointer(); //cluster finished
171 continue;
172 } else if(diff<-fMatch){ //no match
173 break; //insert new cluster
174 }
175 else { //match found, merge it
176 MergeSeq();
177 return;
178 }
179 }
180
181 InsertSeq(); //start new cluster
182}
183
184void AliL3VHDLClusterFinder::MergeSeq()
185{
186#ifdef DEBUG
6aec85a4 187 fprintf(fdeb,"merged with Mean=%d TC=%d (new Merge=%d) %d %d\n",fSeqs[fPList[fRP]].fMean,fSeqs[fPList[fRP]].fTotalCharge,fSeqs[fPList[fRP]].fMerge+1,fRow,fPad);
8e1f5064 188#endif
189 if(fSeqs[fPList[fRP]].fRow==fSeq.fRow){
190 LOG(AliL3Log::kWarning,"AliL3VHDLClusterFinder::","Memory Check")
191 <<"Sequences can be merged on the same rows only."<<ENDLOG;
192 }
193 if(fSeqs[fPList[fRP]].fLastPad+1!=fSeq.fLastPad){
194 LOG(AliL3Log::kWarning,"AliL3VHDLClusterFinder::","Memory Check")
195 <<"Sequences can be merged on consecutive pads only."<<ENDLOG;
196 }
197
198 fSeqs[fPList[fRP]].fMean=fSeq.fMean; //take the new mean
199 fSeqs[fPList[fRP]].fLastPad=fSeq.fLastPad;
200 fSeqs[fPList[fRP]].fTotalCharge+=fSeq.fTotalCharge;
201 fSeqs[fPList[fRP]].fPad+=fSeq.fPad;
202 fSeqs[fPList[fRP]].fTime+=fSeq.fTime;
203 fSeqs[fPList[fRP]].fPad2+=fSeq.fPad2;
204 fSeqs[fPList[fRP]].fTime2+=fSeq.fTime2;
205 fSeqs[fPList[fRP]].fMerge+=1;
206
207 fPList[fWP]=fPList[fRP];
208 fPList[fRP]=N_clmem; //mark it to be free
209
210 IncRPointer();
211 IncWPointer();
212}
213
214void AliL3VHDLClusterFinder::InsertSeq()
215{
216#ifdef DEBUG
6aec85a4 217 fprintf(fdeb,"inserted %d %d\n",fRow,fPad);
8e1f5064 218#endif
219 NextFreeIndex(); //get next index
220 fSeqs[fFirst]=fSeq; //store data
221 fPList[fWP]=fFirst; //store index
222 IncWPointer();
223}
224
225void AliL3VHDLClusterFinder::OutputMemory()
226{
227 Float_t mtime=0,mpad=0;
228 Float_t mtime2=0,mpad2=0;
229 UInt_t tc,row,mno;
230
231 if(fOP!=fFP){
232 //while(fOP!=fFP){
233
234 UInt_t index=fPList[fOP];
235 fPList[fOP]=N_clmem;
236 //cout << fOP << " - " << fFP << " - " << index << endl;
237 IncPointer(fOP);
238 if(index>=N_clmem) return; //nothing to do
239 //if(index>=N_clmem) continue; //nothing to do
240
241 tc=fSeqs[index].fTotalCharge;
242 mno=fSeqs[index].fMerge;
243 row=fSeqs[index].fRow;
244 if(tc!=0){
245 mtime=(Float_t)(fSeqs[index].fTime)/tc;
246 mtime2=sqrt((Float_t)(fSeqs[index].fTime2)/tc-mtime*mtime);
247 }
248 if(tc!=0){
249 mpad=(Float_t)(fSeqs[index].fPad)/tc;
250 mpad2=sqrt((Float_t)(fSeqs[index].fPad2)/tc-mpad*mpad);
251 }
252
253 if(mno<fMinMerge){ //noise cluster
254#ifdef DEBUG
255 fprintf(fdeb,"OutputMemory %d %d: noise cluster (merge): Row %d Pad %.3f Time %.3f TC %d Merge %d\n",fOP,fFP,row,mpad,mtime,tc,mno);
256#endif
257 FreeSeq(index);
258 return;
259 //continue;
260 }
261
262 if(tc<fThreshold){ //total charge below threshold
263#ifdef DEBUG
264 fprintf(fdeb,"OutputMemory %d %d: noise cluster (threshold): Row %d Pad %.3f Time %.3f TC %d Merge %d\n",fOP,fFP,row,mpad,mtime,tc,mno);
265#endif
266 FreeSeq(index);
267 return;
268 //continue;
269 }
270
271#ifdef DEBUG
272 fprintf(fdeb,"OutputMemory %d: Row %d Pad %.3f Time %.3f TC %d Merge %d\n",fNClusters,row,mpad,mtime,tc,mno);
273#endif
274
275 if(stdout)
276 cout<<"WriteCluster: padrow "<<row<<" pad "<<mpad<< " +- "<<mpad2<<" time "<<mtime<<" +- "<<mtime2<<" charge "<<tc<<endl;
277
278 fNClusters++;
279 FreeSeq(index);
280 }
281}
282
283inline void AliL3VHDLClusterFinder::FreeSeq(UShort_t i)
284{
285 ClearSeq(i);
286 IncPointer(fLast,1,N_clmem);
287}
288
289void AliL3VHDLClusterFinder::Clear()
290{
291 fFirst=0; //first list pointer
292 fLast=0; //last list pointer
293
294 fWP=0; //write pointer
295 fRP=0; //read pointer
296 fEP=0; //end read pointer
297 fFP=0; //end flush pointer
298 fOP=0; //output pointer
299
300 fSeq.fRow=0;
301 fSeq.fLastPad=0;
302 fSeq.fTotalCharge=0;
303 fSeq.fPad=0;
304 fSeq.fTime=0;
305 fSeq.fPad2=0;
306 fSeq.fTime2=0;
307 fSeq.fMean=0;
308 fSeq.fMerge=0;
309
310 for(UInt_t i=0;i<N_mem;i++) fPList[i]=N_clmem;
311 for(UInt_t i=0;i<N_clmem;i++) ClearSeq(i);
312}
313
314void AliL3VHDLClusterFinder::ClearSeq(UShort_t i){
315 fSeqs[i].fRow=0;
316 fSeqs[i].fLastPad=0;
317 fSeqs[i].fTotalCharge=0;
318 fSeqs[i].fPad=0;
319 fSeqs[i].fTime=0;
320 fSeqs[i].fPad2=0;
321 fSeqs[i].fTime2=0;
322 fSeqs[i].fMean=0;
323 fSeqs[i].fMerge=0;
324}
325
326void AliL3VHDLClusterFinder::IncPointer(UShort_t &p, Short_t add, UShort_t N){
327 Short_t pp=p;
328 pp+=add;
329 if(pp>=N) p=UShort_t(pp-N);
330 else if(pp<0) p=UShort_t(pp+N);
331 else p=UShort_t(pp);
332}
333
334inline void AliL3VHDLClusterFinder::IncRPointer(){
335 IncPointer(fRP);
336}
337
338inline void AliL3VHDLClusterFinder::IncWPointer(){
339 IncPointer(fWP);
340
341 if(fWP==fOP){
342 LOG(AliL3Log::kWarning,"AliL3VHDLClusterFinder::IncWPointer","Memory Check")
343 <<"Write pointer overwrites output pointer."<<ENDLOG;
344 }
345}
346
347inline void AliL3VHDLClusterFinder::NextFreeIndex(){
348 IncPointer(fFirst,1,N_clmem);
349 if(fFirst==fLast) {
350 LOG(AliL3Log::kFatal,"AliL3VHDLClusterFinder::GetFreeIndex","Memory Check")
351 <<"No space left in sequence list: "<<fFirst<<"=="<<fLast<<ENDLOG;
352 }
353}
354
355#if 0
356void AliL3ClustFinderNew::WriteClusters(Int_t n_clusters,ClusterData *list)
357{
358 Int_t thisrow,thissector;
359 UInt_t counter = fNClusters;
360
361 for(int j=0; j<n_clusters; j++)
362 {
363 if(!list[j].fFlags) continue; //discard 1 pad clusters
364 if(list[j].fTotalCharge < fThreshold) continue; //noise cluster
365
366 Float_t xyz[3];
367 Float_t fpad =(Float_t)list[j].fPad /(Float_t)list[j].fTotalCharge;
368 Float_t fpad2=fXYErr;
369 Float_t ftime =(Float_t)list[j].fTime /(Float_t)list[j].fTotalCharge;
370 Float_t ftime2=fZErr;
371
372 if(fcalcerr) {
373 fpad2=(Float_t)list[j].fPad2/(Float_t)list[j].fTotalCharge - fpad*fpad;
374 fpad2 = sqrt(fpad2);
375 ftime2=(Float_t)list[j].fTime2/(Float_t)list[j].fTotalCharge - ftime*ftime;
376 ftime2 = sqrt(ftime2);
377 }
378
379 if(fstdout==kTRUE)
380 cout<<"WriteCluster: padrow "<<fCurrentRow<<" pad "<<fpad << "+-"<<fpad2<<" time "<<ftime<<"+-"<<ftime2<<" charge "<<list[j].fTotalCharge<<endl;
381
382 AliL3Transform::Slice2Sector(fCurrentSlice,fCurrentRow,thissector,thisrow);
383 AliL3Transform::Raw2Local(xyz,thissector,thisrow,fpad,ftime);
384 if(xyz[0]==0) LOG(AliL3Log::kError,"AliL3ClustFinder","Cluster Finder")
385 <<AliL3Log::kDec<<"Zero cluster"<<ENDLOG;
386 if(fNClusters >= fMaxNClusters)
387 {
388 LOG(AliL3Log::kError,"AliL3ClustFinder::WriteClusters","Cluster Finder")
389 <<AliL3Log::kDec<<"Too many clusters"<<ENDLOG;
390 return;
391 }
392 fSpacePointData[counter].fCharge = list[j].fTotalCharge;
393 fSpacePointData[counter].fX = xyz[0];
394 fSpacePointData[counter].fY = xyz[1];
395 fSpacePointData[counter].fZ = xyz[2];
396 fSpacePointData[counter].fPadRow = fCurrentRow;
397 fSpacePointData[counter].fXYErr = fpad2;
398 fSpacePointData[counter].fZErr = ftime2;
399 fSpacePointData[counter].fID = counter
400 +((fCurrentSlice&0x7f)<<25)+((fCurrentPatch&0x7)<<22);//Uli
401
402 fNClusters++;
403 counter++;
404 }
405}
406#endif
407