Changes done to make the Cluser Finder calculate the errors in Pad and Time direction...
[u/mrichter/AliRoot.git] / HLT / src / AliL3ClustFinder.cxx
1 //*-- Author : Anders Strand Vestbo <mailto:vestbo@fi.uib.no>
2 // @(#) 16.12.2000 
3 //*-- Copyright &copy ASV 
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 //
14 // AliL3ClusterFinder
15 //
16
17
18 ClassImp(AliL3ClustFinder)
19
20 AliL3ClustFinder::AliL3ClustFinder()
21 {
22  
23   fDeconvTime=false;
24   fDeconvPad=false;
25   fXYErr = 0;
26   fZErr = 0;
27 }
28
29 AliL3ClustFinder::AliL3ClustFinder(AliL3Transform *transform)
30 {
31
32   fTransform = transform;
33
34   fThreshold = 20;
35   fDeconvTime = true;
36   fDeconvPad = true;
37 }
38
39
40 AliL3ClustFinder::~AliL3ClustFinder()
41 {
42   //destructor 
43 }
44
45
46 void AliL3ClustFinder::SetTransformer( AliL3Transform *transform )
47 {
48   fTransform = transform;
49 }
50
51
52 void 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
63 void 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
71 void AliL3ClustFinder::SetOutputArray(AliL3SpacePointData *pt)
72 {
73   
74   fSpacePointData = pt;
75 }
76
77 void AliL3ClustFinder::Read(UInt_t ndigits,AliL3DigitRowData *ptr)
78 {
79   fNDigitRowData = ndigits;
80   fDigitRowData = ptr;
81
82 }
83
84 void 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
105 void 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
402 void 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
412 void 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
430 void 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       
459       fTransform->Slice2Sector(fCurrentSlice,fCurrentRow,thissector,thisrow);
460       //if(fCurrentRow > 54) {thisrow = fCurrentRow-55; thissector = fCurrentSlice+36;}
461       //else {thisrow = fCurrentRow; thissector = fCurrentSlice;}
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
478                   +((fCurrentSlice&0x7f)<<25)+((fCurrentPatch&0x7)<<22);
479       
480       fNClusters++;
481       counter++;
482
483     }
484
485   
486 }