8e1f5064 |
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 | |