Using TMath::Abs instead of fabs
[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 "AliL3StandardIncludes.h"
7
8 #include "AliL3VHDLClusterFinder.h"
9
10 #if GCCVERSION == 3
11 using namespace std;
12 #endif
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
29 ClassImp(AliL3VHDLClusterFinder)
30
31 AliL3VHDLClusterFinder::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
52 AliL3VHDLClusterFinder::~AliL3VHDLClusterFinder()
53 {
54 #ifdef DEBUG
55   fclose(fdeb);
56 #endif
57 }
58
59 void 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];
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 #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     fTC=0,fMT=0,fST=0;
88
89     //calculate sequence values 
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;
124     }
125
126     MakeSequence();
127     ProcessSequence();
128     OutputMemory();      
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
141 void 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
162 void 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
179 void 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
188 void 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
198 void AliL3VHDLClusterFinder::CompareSeq()
199 {
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       
208     } else { //match found, merge it
209       MergeSeq(); 
210       return;
211     }
212   }
213
214   InsertSeq(); //start new cluster
215 }
216
217 void AliL3VHDLClusterFinder::MergeSeq()
218 {
219 #ifdef DEBUG
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);
221 #endif
222   if(fSeqs[fPList[fRP]].fRow!=fSeq.fRow){
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   }
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;
240
241     fSeqs[fPList[fRP]].fLastCharge = fSeq.fTotalCharge;
242   }
243   
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
260 void AliL3VHDLClusterFinder::InsertSeq()
261 {
262 #ifdef DEBUG
263   fprintf(fdeb,"inserted %d %d\n",fRow,fPad);
264 #endif
265   NextFreeIndex();    //get next index
266   fSeqs[fFirst]=fSeq; //store data
267   fPList[fWP]=fFirst; //store index
268   IncWPointer();
269 }
270
271 void 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
329 inline void AliL3VHDLClusterFinder::FreeSeq(UShort_t i)
330 {
331   ClearSeq(i);
332   IncPointer(fLast,1,N_clmem);
333 }
334
335 void 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
360 void 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
372 void 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
380 inline void AliL3VHDLClusterFinder::IncRPointer(){
381   IncPointer(fRP);
382 }
383
384 inline 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
393 inline 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
402 void 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