]>
Commit | Line | Data |
---|---|---|
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 |