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