Changes done to make the Cluser Finder calculate the errors in Pad and Time direction...
[u/mrichter/AliRoot.git] / HLT / src / AliL3ClustFinder.cxx
CommitLineData
b661165c 1//*-- Author : Anders Strand Vestbo <mailto:vestbo@fi.uib.no>
2// @(#) 16.12.2000
3//*-- Copyright &copy ASV
108615fc 4
5#include <iostream.h>
6
7#include "AliL3Logging.h"
8#include "AliL3ClustFinder.h"
9#include "AliL3Transform.h"
10#include "AliL3DigitData.h"
11#include "AliL3SpacePointData.h"
12
13//
c3dd27a3 14// AliL3ClusterFinder
108615fc 15//
16
b661165c 17
108615fc 18ClassImp(AliL3ClustFinder)
19
20AliL3ClustFinder::AliL3ClustFinder()
21{
22
23 fDeconvTime=false;
24 fDeconvPad=false;
25 fXYErr = 0;
26 fZErr = 0;
27}
28
29AliL3ClustFinder::AliL3ClustFinder(AliL3Transform *transform)
30{
31
32 fTransform = transform;
33
34 fThreshold = 20;
35 fDeconvTime = true;
36 fDeconvPad = true;
37}
38
39
40AliL3ClustFinder::~AliL3ClustFinder()
41{
42 //destructor
43}
44
45
c3dd27a3 46void AliL3ClustFinder::SetTransformer( AliL3Transform *transform )
47{
48 fTransform = transform;
49}
108615fc 50
51
52void AliL3ClustFinder::InitSlice(Int_t slice,Int_t patch,Int_t firstrow, Int_t lastrow,Int_t nmaxpoints)
53{
54 fNClusters = 0;
55 fMaxNClusters = nmaxpoints;
56 fCurrentSlice = slice;
57 fCurrentPatch = patch;
58 fFirstRow = firstrow;
59 fLastRow = lastrow;
60}
61
62
63void AliL3ClustFinder::InitSlice(Int_t slice,Int_t patch,Int_t nmaxpoints)
64{
65 fNClusters = 0;
66 fMaxNClusters = nmaxpoints;
67 fCurrentSlice = slice;
68 fCurrentPatch = patch;
69}
70
71void AliL3ClustFinder::SetOutputArray(AliL3SpacePointData *pt)
72{
73
74 fSpacePointData = pt;
75}
76
77void AliL3ClustFinder::Read(UInt_t ndigits,AliL3DigitRowData *ptr)
78{
79 fNDigitRowData = ndigits;
80 fDigitRowData = ptr;
81
82}
83
84void AliL3ClustFinder::ProcessDigits()
85{
86 //Loop over rows, and call processrow
87
88
89 AliL3DigitRowData *tempPt = (AliL3DigitRowData*)fDigitRowData;
90
91 for(Int_t i=fFirstRow; i<=fLastRow; i++)
92 {
93 fCurrentRow = i;
94 ProcessRow(tempPt);
95 Byte_t *tmp = (Byte_t*)tempPt;
96 Int_t size = sizeof(AliL3DigitRowData) + tempPt->fNDigit*sizeof(AliL3DigitData);
97 tmp += size;
98 tempPt = (AliL3DigitRowData*)tmp;
99 }
100
101 LOG(AliL3Log::kInformational,"AliL3ClustFinder","Cluster Finder")
102 <<AliL3Log::kDec<<"ClusterFinder found "<<fNClusters<<" clusters"<<ENDLOG;
103}
104
105void AliL3ClustFinder::ProcessRow(AliL3DigitRowData* tempPt)
106{
107
108 Int_t cl_found,trackID[3];
109 //UInt_t rows;
110
111 //Reserve array for found clusters:
112 resx aliresx[5000];
113 bzero((char *)aliresx,5000*sizeof(resx));
114
115 //local results
116 resx *r;
117// UInt_t pres1[2500],pres2[2500];
118 resx *pres1[2500];
119 resx *pres2[2500];
120 resx rr_local;
121 resx *rr = &rr_local;
122
123 //Set variables for each row:
124 Int_t pres_cou1,pres_cou2;
125// UInt_t *r1,*r2;
126 resx **r1,**r2;
127 pres_cou1=pres_cou2=0;
128 r1=pres2;
129 r2=pres1;
130 r=aliresx;
131
132 //variables for each pad:
133 UInt_t start;
134 Int_t start_new=-1;
135 Int_t go_on;
136 Int_t cl_counter=0,rows;
137 Int_t prev_start;
138
139 UInt_t last_pad = 123456789;
140
141 //Loop over digits in this row:
142 for(start=0; start<tempPt->fNDigit; start++)
143 {
144 AliL3DigitData *bins = tempPt->fDigitData;
145
146 if(bins[start].fPad!=last_pad)
147 {//new pad
148
149 last_pad=bins[start].fPad;
150
151 if(r2==pres2)
152 {
153 r1=pres2;
154 r2=pres1;
155 }
156 else
157 {
158 r1=pres1;
159 r2=pres2;
160 }
161 pres_cou1=pres_cou2;
162 pres_cou2=0;
163
164 start_new=-1;
165
166 cl_counter=0;
167 prev_start=0;
168
169
170 }//Looping over pads
171
172 UInt_t av,charge;
173 UInt_t mean;
174 UInt_t flags;
175 Int_t last_falling=0;
176// UInt_t *ri;
177 resx **ri;
178 UInt_t tmp_charge;
179
180 //Set flag:
181 flags=0;
182
183 cl_counter++;
184 if(cl_counter>MAX_C) {LOG(AliL3Log::kWarning,"AliL3ClustFinder::ProcessRow","Cluster Finder")
185 <<AliL3Log::kDec<<"Too many clusters"<<ENDLOG; return;} //too many clusters
186
187
188 if(fDeconvTime)
189 {
190
191 //this is a goto label:
192 redo: ;
193 //divide sequenz in two in case there is a change in slope:
194 if(start_new>-1)
195 {
196 start=start_new;
197 start_new=-1;
198 }
199 }
200 //average and charge
201 av=charge=0;
202
203 //ok, create a new cluster:
204 go_on=1;
205
206 if(fDeconvTime)
207 {
208
209 last_falling=flags=0;
210 }
211
212 if(fDeconvPad)
213 {
214 flags=0;
215 }
216
217
218 //block to introduce new variables:
219 {
220 UInt_t last_a;
221 last_a=0;
222
223 //Loop over this sequenz:
224 while(1)
225 {
226 UInt_t aa;
227 //Int_t start_temp=start;
228
229 //get the adc-value:
230 aa=bins[start].fCharge;
231
232
233 if(fDeconvTime)
234 {
235 //Check if the last pixel in this sequenz is bigger or smaller than this:
236 if(aa>last_a)
237 {
238 if(last_falling)
239 {
240 start_new=start;
241 break;
242 }
243 }
244 else last_falling=1;
245 last_a=aa;
246
247 }
248 //sum of total charge of this cluster on this pad
249 charge +=aa;
250 //this one is needed to determine the mean over all pads:
251 //av += start * (aa);
252 av+=bins[start].fTime*(aa);
253
254 if((bins[start+1].fTime)!=(bins[start].fTime+1) ||
255 ((bins[start+1].fPad)!=(bins[start].fPad))) break;
256
257 start++;
258 }//Loop over sequenz
259 }//block of new variables
260
261 if(fDeconvTime)
262 {
263 if(start_new>0) flags = FLAG_DOUBLE_T;
264 }
265 //calculate mean in time direction for this sequenz:
266 if(charge)
267 {
268 //mean=cog for this sequenz:
269 mean=av/charge;
270 }
271 else
272 {
273 charge=1;
274 mean=1;
275 }
276 //get the pointer to results on previous pad(s)
277 ri=r1;
278
279 //Center of gravity in pad direction
280 //tmp_charge=charge*bins[start-1].fPad;
281 tmp_charge=charge*bins[start].fPad;
282
283 //compare with prevously stored results (previous pads)
284 for(int k=0; k<pres_cou1; k++)
285 {
286 //variable to store difference:
287 Int_t v;
288
289 //structure that contains means (+more) of last pads(s)
290 resx *rr_tmp;
291
292 //go get it
293 rr_tmp = (resx *) *ri;
294
295 //difference between mean and mean on previous pads
296 v=mean-rr_tmp->mean;
297
298
299 if(v<-PARAM1) break;
300
301 //goto next mean on previous pad
302 ri++;
303
304 //is there a matching cluster on previous pad:
305 if(v<=PARAM1)
306 {
307 //ok, we got a new one, so we will add the new and the old
308 rr=&rr_local;
309
310 //Add this function somewhere
311 memcpy7((UInt_t *)rr,(UInt_t *)rr_tmp);
312
313 /* in case we found another sequenz-section on this pad
314 don't loop over all sequenzes on last pad but start
315 where we found a matching candidate for this sequenz-
316 section, adjust upper limit (=pres_cou1)*/
317 r1=ri;
318 pres_cou1-=1+k;
319
320 if(fDeconvPad)
321 {
322 if(charge>rr->scharge)
323 {
324 if(rr->falling)
325 {
326 //previous state was falling
327 flags |= FLAG_DOUBLE_PAD;
328 rr_tmp->flags |= flags;
329 //and create a new one
330 break;
331 }
332 }
333 else
334 {
335 rr->falling = 1;
336 }
337
338 }
339 //don't create a new one
340 go_on=0;
341
342 //store the pointer to the matched in "2"
343// r2[pres_cou2++] = (UInt_t) rr_tmp;
344 r2[pres_cou2++] = rr_tmp;
345
346 //calculate and fill the new means, charge etc. in rr_tmp
347 //charge on this pad to determine whether we cut in pad direction
348 rr->scharge = charge;
349 //unset FLAG_ONEPAD since we have at least 2 pads
350 rr->flags &= (~FLAG_ONEPAD);
351 //summ the charges for de/dx
352 rr->charge += charge;
353 //calculate mean in pad direction
354 rr->pad += tmp_charge;
355 //calculate mean in time direction
356 rr->t+=av;
357 //store the mean of this sequenz on this pad
358 rr->mean = mean;
359
360 //Copy it back to storage
361 memcpy7((UInt_t *)rr_tmp,(UInt_t *)rr);
362
363 //we found a matching one so stop looping over pads on previous pad
364 break;
365
366 }//matching cluster on previous pads
367 }//loop: means on previous pad
368
369 //here we create from scratch new clusters
370 if(go_on)
371 {
372 //store this one beacuse it is the first
373// r2[pres_cou2++]= (UInt_t )r;
374 r2[pres_cou2++]= r;
375
376 mstore((UInt_t *)r,av,tmp_charge,charge,flags|FLAG_ONEPAD,mean,trackID);
377 r++;
378 }
379
380 if(fDeconvTime)
381 {
382
383 //maybe we still have to do something
384 if(start_new>=0) goto redo;
385 }
386
387 }//end of loop over digits
388
389
390 //number of clusters found
391 cl_found = r-aliresx;
392
393 //tot_cl_found+=cl_found;
394
395 //there was something on this padrow
396 rows++;
397
398 WriteClusters(cl_found,aliresx);
399
400}
401
402void AliL3ClustFinder::memcpy7(UInt_t *dst, UInt_t *src)
403{
404 Int_t i;
405 for(i=0; i<7; i++)
406 {
407 *dst++ = *src++;
408 }
409 return;
410}
411
412void AliL3ClustFinder::mstore(UInt_t *r,UInt_t av,UInt_t pad,UInt_t ch,UInt_t flags,UInt_t mean,Int_t *trackID)
413{
414 resx *rr = (resx *) r;
415
416 rr->mean = mean;
417 rr->charge = ch;
418 rr->flags =flags;
419 rr->scharge = ch;
420 rr->pad = pad;
421 rr->t = av;
422 rr->falling=0;
423
424 rr->trackID[0]=trackID[0];
425 rr->trackID[1]=trackID[1];
426 rr->trackID[2]=trackID[2];
427 return;
428}
429
430void AliL3ClustFinder::WriteClusters(Int_t ncl,resx *r)
431{
432/*
433 if(fNClusters >= fMaxNClusters)
434 {
435 LOG(AliL3Log::kError,"AliL3ClustFinder::WriteClusters","Cluster Finder")
436 <<AliL3Log::kDec<<"Too many clusters"<<ENDLOG;
437 return;
438 }
439*/
440 Int_t thisrow,thissector;
441 UInt_t counter = fNClusters;
442
443 if(!fXYErr || !fZErr)
444 {
445 LOG(AliL3Log::kError,"AliL3ClustFinder::WriteClusters","Cluster Errors")
446 <<"No Errors"<<ENDLOG;
447 return;
448 }
449
450 for(int j=0; j<ncl; j++)
451 {
452 if(r[j].flags & FLAG_ONEPAD) continue; //discard 1 pad clusters
453 if(r[j].charge < fThreshold) continue; //noise cluster
454 Float_t xyz[3];
455
456 Float_t fpad=(Float_t)r[j].pad/(Float_t)r[j].charge;
457 Float_t ftime=(Float_t)r[j].t/(Float_t)r[j].charge;
458
ae97a0b9 459 fTransform->Slice2Sector(fCurrentSlice,fCurrentRow,thissector,thisrow);
460 //if(fCurrentRow > 54) {thisrow = fCurrentRow-55; thissector = fCurrentSlice+36;}
461 //else {thisrow = fCurrentRow; thissector = fCurrentSlice;}
108615fc 462 fTransform->Raw2Local(xyz,thissector,thisrow,fpad,ftime);
463 if(xyz[0]==0) LOG(AliL3Log::kError,"AliL3ClustFinder","Cluster Finder")
464 <<AliL3Log::kDec<<"Zero cluster"<<ENDLOG;
465 if(fNClusters >= fMaxNClusters)
466 {
467 LOG(AliL3Log::kError,"AliL3ClustFinder::WriteClusters","Cluster Finder")
468 <<AliL3Log::kDec<<"Too many clusters"<<ENDLOG;
469 return;
470 }
471 fSpacePointData[counter].fX = xyz[0];
472 fSpacePointData[counter].fY = xyz[1];
473 fSpacePointData[counter].fZ = xyz[2];
474 fSpacePointData[counter].fPadRow = fCurrentRow;
475 fSpacePointData[counter].fXYErr = fXYErr;
476 fSpacePointData[counter].fZErr = fZErr;
477 fSpacePointData[counter].fID = counter
34c9ed82 478 +((fCurrentSlice&0x7f)<<25)+((fCurrentPatch&0x7)<<22);
108615fc 479
480 fNClusters++;
481 counter++;
482
483 }
484
485
486}