Added AliL3Stopwatch.
[u/mrichter/AliRoot.git] / HLT / misc / AliL3VHDLClusterFinder.cxx
1 // $Id$
2
3 // Author: Constantin Loizides <mailto:loizides@ikf.physik.uni-frankfurt.de>
4 //*-- Copyright & copy CL
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
25 //</pre> 
26 */
27
28 ClassImp(AliL3VHDLClusterFinder)
29
30 AliL3VHDLClusterFinder::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
51 AliL3VHDLClusterFinder::~AliL3VHDLClusterFinder()
52 {
53 #ifdef DEBUG
54   fclose(fdeb);
55 #endif
56 }
57
58 void 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;
105     //if(tc!=0) fSeq.fMean=mt/tc;
106     if(tc!=0) fSeq.fMean=rtime-i/2;    
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
127 void 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
144 void 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
153 void 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
163 void 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
184 void AliL3VHDLClusterFinder::MergeSeq()
185 {
186 #ifdef DEBUG
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);
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
214 void AliL3VHDLClusterFinder::InsertSeq()
215 {
216 #ifdef DEBUG
217   fprintf(fdeb,"inserted %d %d\n",fRow,fPad);
218 #endif
219   NextFreeIndex();    //get next index
220   fSeqs[fFirst]=fSeq; //store data
221   fPList[fWP]=fFirst; //store index
222   IncWPointer();
223 }
224
225 void 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
283 inline void AliL3VHDLClusterFinder::FreeSeq(UShort_t i)
284 {
285   ClearSeq(i);
286   IncPointer(fLast,1,N_clmem);
287 }
288
289 void 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
314 void 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
326 void 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
334 inline void AliL3VHDLClusterFinder::IncRPointer(){
335   IncPointer(fRP);
336 }
337
338 inline 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
347 inline 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
356 void 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