Created cluster finder class that simulates the VHDL cluster finder on Altro data.
[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 ClassImp(AliL3VHDLClusterFinder)
27
28 AliL3VHDLClusterFinder::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
49 AliL3VHDLClusterFinder::~AliL3VHDLClusterFinder()
50 {
51 #ifdef DEBUG
52   fclose(fdeb);
53 #endif
54 }
55
56 void 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
124 void 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
141 void 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
150 void 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
160 void 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
181 void 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
211 void 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
222 void 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
280 inline void AliL3VHDLClusterFinder::FreeSeq(UShort_t i)
281 {
282   ClearSeq(i);
283   IncPointer(fLast,1,N_clmem);
284 }
285
286 void 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
311 void 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
323 void 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
331 inline void AliL3VHDLClusterFinder::IncRPointer(){
332   IncPointer(fRP);
333 }
334
335 inline 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
344 inline 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
353 void 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