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