]>
Commit | Line | Data |
---|---|---|
1 | // @(#) $Id$ | |
2 | ||
3 | // Author: Constantin Loizides <mailto:loizides@ikf.uni-frankfurt.de> | |
4 | //*-- Copyright & copy ALICE HLT Group | |
5 | /** \class AliL3VHDLClusterFinder | |
6 | <pre> | |
7 | //____________________________________________________ | |
8 | // AliL3VHDLClusterFinder | |
9 | // | |
10 | // The current VHDL cluster finder for HLT | |
11 | // Based on STAR L3 | |
12 | // | |
13 | // Most important parameters: | |
14 | // fThreshold - threshold for noise clusters | |
15 | // fMatch - length in time for overlapping sequences | |
16 | </pre> | |
17 | */ | |
18 | ||
19 | #include "AliL3StandardIncludes.h" | |
20 | #include "AliL3RootTypes.h" | |
21 | #include "AliL3Logging.h" | |
22 | #include "AliL3AltroMemHandler.h" | |
23 | ||
24 | //#define VHDLDEBUG | |
25 | #include "AliL3VHDLClusterFinder.h" | |
26 | ||
27 | #if __GNUC__ >= 3 | |
28 | using namespace std; | |
29 | #endif | |
30 | ||
31 | ||
32 | ClassImp(AliL3VHDLClusterFinder) | |
33 | ||
34 | AliL3VHDLClusterFinder::AliL3VHDLClusterFinder() | |
35 | { | |
36 | // default constructor | |
37 | fMatch = 4; | |
38 | fThreshold = 10; | |
39 | fMinMerge = 1; | |
40 | fNClusters=0; | |
41 | ||
42 | fXYErr = 0.2; | |
43 | fZErr = 0.3; | |
44 | fDeconvPad = kTRUE; | |
45 | fDeconvTime = kTRUE; | |
46 | fstdout = kFALSE; | |
47 | fcalcerr = kTRUE; | |
48 | ||
49 | Clear(); | |
50 | #ifdef VHDLDEBUG | |
51 | fdeb=fopen("vhdlclusterfinder.debug","w"); | |
52 | //fdeb=stderr; | |
53 | #endif | |
54 | } | |
55 | ||
56 | AliL3VHDLClusterFinder::~AliL3VHDLClusterFinder() | |
57 | { | |
58 | // destructor | |
59 | #ifdef VHDLDEBUG | |
60 | fclose(fdeb); | |
61 | #endif | |
62 | } | |
63 | ||
64 | void AliL3VHDLClusterFinder::ProcessDigits() | |
65 | { | |
66 | //Loop over data like the analyzer of the VHDL code | |
67 | const UChar_t kn=255; | |
68 | UShort_t rrow=0,rtime=0; | |
69 | UChar_t rpad=0,i=kn; | |
70 | UShort_t *charges=new UShort_t[kn]; | |
71 | ||
72 | fNClusters=0; | |
73 | fRow=0; | |
74 | fPad=0; | |
75 | Clear(); | |
76 | ||
77 | //loop over input data | |
78 | while(fAltromem.ReadSequence(rrow,rpad,rtime,i,&charges)){ | |
79 | #if 0 | |
80 | cout << "Padrow " << (int)rrow << " pad " << (int)rpad << " time " <<(int)rtime << " charges "; | |
81 | for(UChar_t ii=0;ii<i;ii++) cout << (int)charges[ii] << " "; | |
82 | cout << endl; | |
83 | #endif | |
84 | #ifdef VHDLDEBUG | |
85 | fprintf(fdeb,"ProcessDigits: Input Data: %d %d %d charges:",(int)rrow,(int)rpad,(int)rtime); | |
86 | for(UChar_t ii=0;ii<i;ii++) fprintf(fdeb," %d",(int)charges[ii]); | |
87 | fprintf(fdeb,"\n"); | |
88 | #endif | |
89 | ||
90 | fNRow=rrow; | |
91 | fNPad=rpad; | |
92 | fTC=0,fMT=0,fST=0; | |
93 | ||
94 | //calculate sequence values | |
95 | if(fDeconvTime){ | |
96 | UChar_t ii=0; | |
97 | Int_t sl=0; | |
98 | Int_t charge=0,lcharge=0; | |
99 | Bool_t falling=kFALSE; | |
100 | while(ii<i){ | |
101 | charge=charges[ii]; | |
102 | if((falling)&&(charge>lcharge)){ | |
103 | fSM=rtime-sl/2; | |
104 | MakeSequence(); | |
105 | ProcessSequence(); | |
106 | OutputMemory(); | |
107 | ||
108 | rtime=rtime-sl; | |
109 | falling=kFALSE; | |
110 | sl=0; | |
111 | fTC=0,fMT=0,fST=0; | |
112 | } else if(charge<lcharge) falling=kTRUE; | |
113 | ||
114 | fTC+=charge; | |
115 | fMT+=(rtime-sl)*charge; | |
116 | fST+=(rtime-sl)*(rtime-sl)*charge; | |
117 | sl++; | |
118 | ii++; | |
119 | lcharge=charge; | |
120 | } //end loop over sequence | |
121 | fSM=rtime-sl/2; | |
122 | } else { /* no deconvolution */ | |
123 | for(UChar_t ii=0;ii<i;ii++){ | |
124 | fTC+=charges[ii]; | |
125 | fMT+=(rtime-ii)*charges[ii]; | |
126 | fST+=(rtime-ii)*(rtime-ii)*charges[ii]; | |
127 | } | |
128 | fSM=rtime-i/2; | |
129 | } | |
130 | ||
131 | MakeSequence(); | |
132 | ProcessSequence(); | |
133 | OutputMemory(); | |
134 | ||
135 | #ifdef VHDLDEBUG | |
136 | fflush(fdeb); | |
137 | #endif | |
138 | i=kn; //store size of charges array | |
139 | } //loop over data | |
140 | ||
141 | //flush everything left | |
142 | FlushMemory(); | |
143 | while(fOP!=fFP) OutputMemory(); | |
144 | } | |
145 | ||
146 | void AliL3VHDLClusterFinder::MakeSequence(){ | |
147 | // makes the sequence | |
148 | if(!fTC) return; | |
149 | ||
150 | Int_t mp=fNPad*fTC; | |
151 | Int_t sp=fNPad*fNPad*fTC; | |
152 | ||
153 | fSeq.fTotalCharge=fTC; | |
154 | fSeq.fPad=mp; | |
155 | fSeq.fTime=fMT; | |
156 | fSeq.fPad2=sp; | |
157 | fSeq.fTime2=fST; | |
158 | fSeq.fMean=0; | |
159 | fSeq.fMean=fSM; | |
160 | fSeq.fMerge=0; | |
161 | fSeq.fRow=fNRow; | |
162 | fSeq.fLastPad=fNPad; | |
163 | fSeq.fChargeFalling=0; | |
164 | if(fDeconvPad) fSeq.fLastCharge=fTC; | |
165 | else fSeq.fLastCharge=0; | |
166 | } | |
167 | ||
168 | void AliL3VHDLClusterFinder::ProcessSequence() | |
169 | { | |
170 | // processes the sequence | |
171 | if(fNRow!=fRow) FlushMemory(); | |
172 | else if(fNPad==fPad+1) PrepareMemory(); | |
173 | else if(fNPad!=fPad) FlushMemory(); | |
174 | ||
175 | //store new row and pad values | |
176 | fRow=fNRow; | |
177 | fPad=fNPad; | |
178 | ||
179 | #ifdef VHDLDEBUG | |
180 | fprintf(fdeb,"ProcessSequence: Mean=%d TC=%d ",fSeq.fMean,fSeq.fTotalCharge); | |
181 | #endif | |
182 | ||
183 | CompareSeq(); //merge or insert | |
184 | } | |
185 | ||
186 | void AliL3VHDLClusterFinder::PrepareMemory() | |
187 | { | |
188 | // prepares the memory | |
189 | #ifdef VHDLDEBUG | |
190 | fprintf(fdeb,"PrepareMemory %d %d %d\n",fRP,fEP,fWP); | |
191 | #endif | |
192 | fRP=fEP; | |
193 | fEP=fWP; | |
194 | } | |
195 | ||
196 | void AliL3VHDLClusterFinder::FlushMemory() | |
197 | { | |
198 | // flushes the memory | |
199 | #ifdef VHDLDEBUG | |
200 | fprintf(fdeb,"FlushMemory %d %d %d %d\n",fFP,fRP,fEP,fWP); | |
201 | #endif | |
202 | fFP=fWP; | |
203 | fRP=fWP; | |
204 | fEP=fWP; | |
205 | } | |
206 | ||
207 | void AliL3VHDLClusterFinder::CompareSeq() | |
208 | { | |
209 | // compares sequences | |
210 | while(fRP!=fEP){ | |
211 | Int_t diff=fSeqs[fPList[fRP]].fMean-fSeq.fMean; | |
212 | ||
213 | if(diff>fMatch){ //no match | |
214 | IncRPointer(); //cluster finished | |
215 | continue; | |
216 | } else if(diff<-fMatch){ //no match | |
217 | break; //insert new cluster | |
218 | } else { //match found, merge it | |
219 | MergeSeq(); | |
220 | return; | |
221 | } | |
222 | } | |
223 | ||
224 | InsertSeq(); //start new cluster | |
225 | } | |
226 | ||
227 | void AliL3VHDLClusterFinder::MergeSeq() | |
228 | { | |
229 | // merges sequences | |
230 | #ifdef VHDLDEBUG | |
231 | 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); | |
232 | #endif | |
233 | if(fSeqs[fPList[fRP]].fRow!=fSeq.fRow){ | |
234 | LOG(AliL3Log::kWarning,"AliL3VHDLClusterFinder::","Memory Check") | |
235 | <<"Sequences can be merged on the same rows only."<<ENDLOG; | |
236 | } | |
237 | if(fSeqs[fPList[fRP]].fLastPad+1!=fSeq.fLastPad){ | |
238 | LOG(AliL3Log::kWarning,"AliL3VHDLClusterFinder::","Memory Check") | |
239 | <<"Sequences can be merged on consecutive pads only."<<ENDLOG; | |
240 | } | |
241 | ||
242 | if(fDeconvPad){ | |
243 | if(fSeq.fTotalCharge > fSeqs[fPList[fRP]].fLastCharge){ | |
244 | if(fSeqs[fPList[fRP]].fChargeFalling){ //The previous pad was falling | |
245 | IncRPointer(); | |
246 | InsertSeq(); //start a new cluster | |
247 | return; | |
248 | } | |
249 | } | |
250 | else fSeqs[fPList[fRP]].fChargeFalling = 1; | |
251 | ||
252 | fSeqs[fPList[fRP]].fLastCharge = fSeq.fTotalCharge; | |
253 | } | |
254 | ||
255 | fSeqs[fPList[fRP]].fMean=fSeq.fMean; //take the new mean | |
256 | fSeqs[fPList[fRP]].fLastPad=fSeq.fLastPad; | |
257 | fSeqs[fPList[fRP]].fTotalCharge+=fSeq.fTotalCharge; | |
258 | fSeqs[fPList[fRP]].fPad+=fSeq.fPad; | |
259 | fSeqs[fPList[fRP]].fTime+=fSeq.fTime; | |
260 | fSeqs[fPList[fRP]].fPad2+=fSeq.fPad2; | |
261 | fSeqs[fPList[fRP]].fTime2+=fSeq.fTime2; | |
262 | fSeqs[fPList[fRP]].fMerge+=1; | |
263 | ||
264 | fPList[fWP]=fPList[fRP]; | |
265 | fPList[fRP]=N_clmem; //mark it to be free | |
266 | ||
267 | IncRPointer(); | |
268 | IncWPointer(); | |
269 | } | |
270 | ||
271 | void AliL3VHDLClusterFinder::InsertSeq() | |
272 | { | |
273 | // inserts sequence | |
274 | #ifdef VHDLDEBUG | |
275 | fprintf(fdeb,"inserted %d %d\n",fRow,fPad); | |
276 | #endif | |
277 | NextFreeIndex(); //get next index | |
278 | fSeqs[fFirst]=fSeq; //store data | |
279 | fPList[fWP]=fFirst; //store index | |
280 | IncWPointer(); | |
281 | } | |
282 | ||
283 | void AliL3VHDLClusterFinder::OutputMemory() | |
284 | { | |
285 | // output memory? | |
286 | Float_t mtime=0,mpad=0; | |
287 | Float_t mtime2=0,mpad2=0; | |
288 | UInt_t tc,row,mno; | |
289 | ||
290 | if(fOP!=fFP){ | |
291 | //while(fOP!=fFP){ | |
292 | ||
293 | UInt_t index=fPList[fOP]; | |
294 | fPList[fOP]=N_clmem; | |
295 | //cout << fOP << " - " << fFP << " - " << index << endl; | |
296 | IncPointer(fOP); | |
297 | if(index>=N_clmem) return; //nothing to do | |
298 | //if(index>=N_clmem) continue; //nothing to do | |
299 | ||
300 | tc=fSeqs[index].fTotalCharge; | |
301 | mno=fSeqs[index].fMerge; | |
302 | row=fSeqs[index].fRow; | |
303 | if(tc!=0){ | |
304 | mtime=(Float_t)(fSeqs[index].fTime)/tc; | |
305 | mtime2=sqrt((Float_t)(fSeqs[index].fTime2)/tc-mtime*mtime); | |
306 | } | |
307 | if(tc!=0){ | |
308 | mpad=(Float_t)(fSeqs[index].fPad)/tc; | |
309 | mpad2=sqrt((Float_t)(fSeqs[index].fPad2)/tc-mpad*mpad); | |
310 | } | |
311 | ||
312 | if(mno<fMinMerge){ //noise cluster | |
313 | #ifdef VHDLDEBUG | |
314 | 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); | |
315 | #endif | |
316 | FreeSeq(index); | |
317 | return; | |
318 | //continue; | |
319 | } | |
320 | ||
321 | if(tc<fThreshold){ //total charge below threshold | |
322 | #ifdef VHDLDEBUG | |
323 | 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); | |
324 | #endif | |
325 | FreeSeq(index); | |
326 | return; | |
327 | //continue; | |
328 | } | |
329 | ||
330 | #ifdef VHDLDEBUG | |
331 | fprintf(fdeb,"OutputMemory %d: Row %d Pad %.3f Time %.3f TC %d Merge %d\n",fNClusters,row,mpad,mtime,tc,mno); | |
332 | #endif | |
333 | ||
334 | // if(stdout) | |
335 | // cout<<"WriteCluster: padrow "<<row<<" pad "<<mpad<< " +- "<<mpad2<<" time "<<mtime<<" +- "<<mtime2<<" charge "<<tc<<endl; | |
336 | ||
337 | fNClusters++; | |
338 | FreeSeq(index); | |
339 | } | |
340 | } | |
341 | ||
342 | void AliL3VHDLClusterFinder::FreeSeq(UShort_t i) | |
343 | { | |
344 | // frees the sequence | |
345 | ClearSeq(i); | |
346 | IncPointer(fLast,1,N_clmem); | |
347 | } | |
348 | ||
349 | void AliL3VHDLClusterFinder::Clear() | |
350 | { | |
351 | // clears everything | |
352 | fFirst=0; //first list pointer | |
353 | fLast=0; //last list pointer | |
354 | ||
355 | fWP=0; //write pointer | |
356 | fRP=0; //read pointer | |
357 | fEP=0; //end read pointer | |
358 | fFP=0; //end flush pointer | |
359 | fOP=0; //output pointer | |
360 | ||
361 | fSeq.fRow=0; | |
362 | fSeq.fLastPad=0; | |
363 | fSeq.fTotalCharge=0; | |
364 | fSeq.fPad=0; | |
365 | fSeq.fTime=0; | |
366 | fSeq.fPad2=0; | |
367 | fSeq.fTime2=0; | |
368 | fSeq.fMean=0; | |
369 | fSeq.fMerge=0; | |
370 | ||
371 | for(UInt_t i=0;i<N_mem;i++) fPList[i]=N_clmem; | |
372 | for(UInt_t i=0;i<N_clmem;i++) ClearSeq(i); | |
373 | } | |
374 | ||
375 | void AliL3VHDLClusterFinder::ClearSeq(UShort_t i){ | |
376 | // clears a sequence | |
377 | fSeqs[i].fRow=0; | |
378 | fSeqs[i].fLastPad=0; | |
379 | fSeqs[i].fTotalCharge=0; | |
380 | fSeqs[i].fPad=0; | |
381 | fSeqs[i].fTime=0; | |
382 | fSeqs[i].fPad2=0; | |
383 | fSeqs[i].fTime2=0; | |
384 | fSeqs[i].fMean=0; | |
385 | fSeqs[i].fMerge=0; | |
386 | } | |
387 | ||
388 | void AliL3VHDLClusterFinder::IncPointer(UShort_t &p, Short_t add, UShort_t N){ | |
389 | // increments pointer by 'add' in a circular buffer | |
390 | Short_t pp=p; | |
391 | pp+=add; | |
392 | if(pp>=N) p=UShort_t(pp-N); | |
393 | else if(pp<0) p=UShort_t(pp+N); | |
394 | else p=UShort_t(pp); | |
395 | } | |
396 | ||
397 | void AliL3VHDLClusterFinder::IncRPointer(){ | |
398 | // increments pointer fRP | |
399 | IncPointer(fRP); | |
400 | } | |
401 | ||
402 | void AliL3VHDLClusterFinder::IncWPointer(){ | |
403 | // increments pointer fWP | |
404 | IncPointer(fWP); | |
405 | ||
406 | if(fWP==fOP){ | |
407 | LOG(AliL3Log::kWarning,"AliL3VHDLClusterFinder::IncWPointer","Memory Check") | |
408 | <<"Write pointer overwrites output pointer."<<ENDLOG; | |
409 | } | |
410 | } | |
411 | ||
412 | void AliL3VHDLClusterFinder::NextFreeIndex(){ | |
413 | // finds next free index | |
414 | IncPointer(fFirst,1,N_clmem); | |
415 | if(fFirst==fLast) { | |
416 | LOG(AliL3Log::kFatal,"AliL3VHDLClusterFinder::GetFreeIndex","Memory Check") | |
417 | <<"No space left in sequence list: "<<fFirst<<"=="<<fLast<<ENDLOG; | |
418 | } | |
419 | } | |
420 | ||
421 | #if 0 | |
422 | void AliL3ClustFinderNew::WriteClusters(Int_t n_clusters,ClusterData *list) | |
423 | { | |
424 | // writes clusters | |
425 | Int_t thisrow,thissector; | |
426 | UInt_t counter = fNClusters; | |
427 | ||
428 | for(int j=0; j<n_clusters; j++) | |
429 | { | |
430 | if(!list[j].fFlags) continue; //discard 1 pad clusters | |
431 | if(list[j].fTotalCharge < fThreshold) continue; //noise cluster | |
432 | ||
433 | Float_t xyz[3]; | |
434 | Float_t fpad =(Float_t)list[j].fPad /(Float_t)list[j].fTotalCharge; | |
435 | Float_t fpad2=fXYErr; | |
436 | Float_t ftime =(Float_t)list[j].fTime /(Float_t)list[j].fTotalCharge; | |
437 | Float_t ftime2=fZErr; | |
438 | ||
439 | if(fcalcerr) { | |
440 | fpad2=(Float_t)list[j].fPad2/(Float_t)list[j].fTotalCharge - fpad*fpad; | |
441 | fpad2 = sqrt(fpad2); | |
442 | ftime2=(Float_t)list[j].fTime2/(Float_t)list[j].fTotalCharge - ftime*ftime; | |
443 | ftime2 = sqrt(ftime2); | |
444 | } | |
445 | ||
446 | if(fstdout==kTRUE) | |
447 | cout<<"WriteCluster: padrow "<<fCurrentRow<<" pad "<<fpad << "+-"<<fpad2<<" time "<<ftime<<"+-"<<ftime2<<" charge "<<list[j].fTotalCharge<<endl; | |
448 | ||
449 | AliL3Transform::Slice2Sector(fCurrentSlice,fCurrentRow,thissector,thisrow); | |
450 | AliL3Transform::Raw2Local(xyz,thissector,thisrow,fpad,ftime); | |
451 | if(xyz[0]==0) LOG(AliL3Log::kError,"AliL3ClustFinder","Cluster Finder") | |
452 | <<AliL3Log::kDec<<"Zero cluster"<<ENDLOG; | |
453 | if(fNClusters >= fMaxNClusters) | |
454 | { | |
455 | LOG(AliL3Log::kError,"AliL3ClustFinder::WriteClusters","Cluster Finder") | |
456 | <<AliL3Log::kDec<<"Too many clusters"<<ENDLOG; | |
457 | return; | |
458 | } | |
459 | fSpacePointData[counter].fCharge = list[j].fTotalCharge; | |
460 | fSpacePointData[counter].fX = xyz[0]; | |
461 | fSpacePointData[counter].fY = xyz[1]; | |
462 | fSpacePointData[counter].fZ = xyz[2]; | |
463 | fSpacePointData[counter].fPadRow = fCurrentRow; | |
464 | fSpacePointData[counter].fXYErr = fpad2; | |
465 | fSpacePointData[counter].fZErr = ftime2; | |
466 | fSpacePointData[counter].fID = counter | |
467 | +((fCurrentSlice&0x7f)<<25)+((fCurrentPatch&0x7)<<22);//Uli | |
468 | ||
469 | fNClusters++; | |
470 | counter++; | |
471 | } | |
472 | } | |
473 | #endif | |
474 |