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