b94af6883be6e79b9a1f1d610d43db039a6f2d2d
[u/mrichter/AliRoot.git] / HLT / src / AliL3ClustFinderNew.cxx
1 //$Id$
2
3 // Author: Anders Vestbo <mailto:vestbo@fi.uib.no>
4 //*-- Copyright &copy ASV 
5
6 #include <iostream.h>
7 #include "AliL3Logging.h"
8 #include "AliL3ClustFinderNew.h"
9 #include "AliL3DigitData.h"
10 #include "AliL3Transform.h"
11 #include "AliL3SpacePointData.h"
12
13 //_____________________________________________________________
14 // AliL3ClustFinderNew
15 //
16 // The current cluster finder for HLT
17 // Based on STAR L3
18
19 ClassImp(AliL3ClustFinderNew)
20
21 AliL3ClustFinderNew::AliL3ClustFinderNew()
22 {
23   fMatch = 4;
24   fThreshold =10;
25   fDeconvPad = kTRUE;
26   fDeconvTime = kTRUE;
27 }
28
29 AliL3ClustFinderNew::AliL3ClustFinderNew(AliL3Transform *transform)
30 {
31   fTransform = transform;
32   fMatch = 4;
33   fThreshold =10;
34   fDeconvTime = kTRUE;
35   fDeconvPad = kTRUE;
36 }
37
38 AliL3ClustFinderNew::~AliL3ClustFinderNew()
39 {
40
41
42 }
43
44 void AliL3ClustFinderNew::InitSlice(Int_t slice,Int_t patch,Int_t firstrow, Int_t lastrow,Int_t nmaxpoints)
45 {
46   fNClusters = 0;
47   fMaxNClusters = nmaxpoints;
48   fCurrentSlice = slice;
49   fCurrentPatch = patch;
50   fFirstRow = firstrow;
51   fLastRow = lastrow;
52   fDeconvTime = kTRUE;
53   fDeconvPad = kTRUE;
54 }
55
56 void AliL3ClustFinderNew::InitSlice(Int_t slice,Int_t patch,Int_t nmaxpoints)
57 {
58   fNClusters = 0;
59   fMaxNClusters = nmaxpoints;
60   fCurrentSlice = slice;
61   fCurrentPatch = patch;
62 }
63
64 void AliL3ClustFinderNew::SetOutputArray(AliL3SpacePointData *pt)
65 {
66   
67   fSpacePointData = pt;
68 }
69
70
71 void AliL3ClustFinderNew::Read(UInt_t ndigits,AliL3DigitRowData *ptr)
72 {
73   fNDigitRowData = ndigits;
74   fDigitRowData = ptr;
75   
76 }
77
78 void AliL3ClustFinderNew::ProcessDigits()
79 {
80   //Loop over rows, and call processrow
81   
82   
83   AliL3DigitRowData *tempPt = (AliL3DigitRowData*)fDigitRowData;
84   
85   for(Int_t i=fFirstRow; i<=fLastRow; i++)
86     {
87       fCurrentRow = i;
88       ProcessRow(tempPt);
89       Byte_t *tmp = (Byte_t*)tempPt;
90       Int_t size = sizeof(AliL3DigitRowData) + tempPt->fNDigit*sizeof(AliL3DigitData);
91       tmp += size;
92       tempPt = (AliL3DigitRowData*)tmp;
93     }
94 }
95
96
97
98 void AliL3ClustFinderNew::ProcessRow(AliL3DigitRowData *tempPt)
99 {
100
101   UInt_t last_pad = 123456789;
102
103   ClusterData *pad1[2500]; //2 lists for internal memory=2pads
104   ClusterData *pad2[2500]; //2 lists for internal memory=2pads
105   ClusterData clusterlist[5000]; //Clusterlist
106
107   ClusterData **currentPt; //List of pointers to the current pad
108   ClusterData **previousPt; //List of pointers to the previous pad
109   currentPt = pad2;
110   previousPt = pad1;
111   UInt_t n_previous=0,n_current=0,n_total=0;
112
113   //Loop over sequences of this row:
114   for(UInt_t bin=0; bin<tempPt->fNDigit; bin++)
115     {
116       
117       AliL3DigitData *data = tempPt->fDigitData;
118       if(data[bin].fPad != last_pad)
119         {
120           //This is a new pad
121           
122           //Switch:
123           if(currentPt == pad2)
124             {
125               currentPt = pad1;
126               previousPt = pad2;
127               
128             }
129           else 
130             {
131               currentPt = pad2;
132               previousPt = pad1;
133             }
134           n_previous = n_current;
135           n_current = 0;
136           if(bin[data].fPad != last_pad+1)
137             {
138               //this happens if there is a pad with no signal.
139               n_previous = n_current = 0;
140             }
141           last_pad = data[bin].fPad;
142         }
143
144       Bool_t new_cluster = kTRUE;
145
146       UInt_t seq_charge=0,seq_average=0;
147       
148       UInt_t last_charge=0,last_was_falling=0;
149       Int_t new_bin=-1;
150       if(fDeconvTime)
151         {
152         redo: //This is a goto.
153           if(new_bin > -1)
154             {
155               bin = new_bin;
156               new_bin = -1;
157             }
158           
159           last_charge=0;
160           last_was_falling = 0;
161         }
162       
163
164       while(1) //Loop over current sequence
165         {
166           if(data[bin].fTime >= fTransform->GetNTimeBins())
167             {
168               LOG(AliL3Log::kFatal,"AliL3ClustFinderNew::ProcessRow","Digits")
169                 <<"Timebin out of range "<<(Int_t)data[bin].fTime<<ENDLOG;
170               break;
171             }
172
173           //Get the current ADC-value
174           UInt_t charge = data[bin].fCharge;
175           
176           if(fDeconvTime)
177             {
178               //Check if the last pixel in the sequence is smaller than this
179               if(charge > last_charge)
180                 {
181                   if(last_was_falling)
182                     {
183                       new_bin = bin;
184                       break;
185                     }
186                 }
187               else last_was_falling = 1; //last pixel was larger than this
188               last_charge = charge;
189             }
190           
191           //Sum the total charge of this sequence
192           seq_charge += charge;
193           seq_average += data[bin].fTime*charge;
194           
195           //Check where to stop:
196           if(data[bin+1].fPad != data[bin].fPad) //new pad
197             break; 
198           if(data[bin+1].fTime != data[bin].fTime+1) //end of sequence
199             break;
200           
201           bin++;
202         }//end loop over sequence
203       
204       //Calculate mean of sequence:
205       Int_t seq_mean=0;
206       if(seq_charge)
207         seq_mean = seq_average/seq_charge;
208       else
209         {
210           LOG(AliL3Log::kFatal,"AliL3ClustFinderNew::ProcessRow","Data")
211             <<"Error in data given to the cluster finder"<<ENDLOG;
212           seq_mean = 1;
213           seq_charge = 1;
214         }
215       
216       //Calculate mean in pad direction:
217       Int_t pad_mean = seq_charge*data[bin].fPad;
218
219       //Compare with results on previous pad:
220       for(UInt_t p=0; p<n_previous; p++)
221         {
222           Int_t difference = seq_mean - previousPt[p]->fMean;
223           if(difference < -fMatch)
224             break;
225
226           if(difference <= fMatch)//There is a match here!!
227             {
228               
229               ClusterData *local = previousPt[p];
230               if(fDeconvPad)
231                 {
232                   if(seq_charge > local->fLastCharge)
233                     {
234                       if(local->fChargeFalling)//The previous pad was falling
235                         {                       
236                           break;//create a new cluster
237                         }                   
238                     }
239                   else
240                     local->fChargeFalling = 1;
241                   local->fLastCharge = seq_charge;
242                 }
243               
244               //Don't create a new cluster, because we found a match
245               new_cluster = kFALSE;
246               
247               //Update cluster on current pad with the matching one:
248               //currentPt[n_current] = previousPt[p];
249               
250               local->fTotalCharge += seq_charge;
251               local->fPad += pad_mean;
252               local->fTime += seq_average;
253               local->fMean = seq_mean;
254               local->fFlags++; //means we have more than one pad 
255               
256               currentPt[n_current] = local;
257               n_current++;
258               
259               break;
260             }//Checking for match at previous pad
261         }//Loop over results on previous pad.
262       
263       if(new_cluster)
264         {
265           //Start a new cluster. Add it to the clusterlist, and update
266           //the list of pointers to clusters in current pad.
267           //current pad will be previous pad on next pad.
268
269           //Add to the clusterlist:
270           ClusterData *tmp = &clusterlist[n_total];
271           tmp->fTotalCharge = seq_charge;
272           tmp->fPad = pad_mean;
273           tmp->fTime = seq_average;
274           tmp->fMean = seq_mean;
275           tmp->fFlags = 0;  //flags for 1 pad clusters
276           if(fDeconvPad)
277             {
278               tmp->fChargeFalling = 0;
279               tmp->fLastCharge = seq_charge;
280             }
281
282           //Update list of pointers to previous pad:
283           currentPt[n_current] = &clusterlist[n_total];
284           n_total++;
285           n_current++;
286         }
287       if(fDeconvTime)
288         if(new_bin >= 0) goto redo;
289       
290     }//Loop over digits on this padrow
291   
292   WriteClusters(n_total,clusterlist);
293
294
295 }
296
297 void AliL3ClustFinderNew::WriteClusters(Int_t n_clusters,ClusterData *list)
298 {
299   Int_t thisrow,thissector;
300   UInt_t counter = fNClusters;
301   
302   for(int j=0; j<n_clusters; j++)
303     {
304       if(!list[j].fFlags) continue; //discard 1 pad clusters
305       if(list[j].fTotalCharge < fThreshold) continue; //noise cluster
306       Float_t xyz[3];
307       
308       Float_t fpad=(Float_t)list[j].fPad/(Float_t)list[j].fTotalCharge;
309       Float_t ftime=(Float_t)list[j].fTime/(Float_t)list[j].fTotalCharge;
310       //printf("padrow %d number of pads %d totalcharge %d\n",fCurrentRow,list[j].fFlags,list[j].fTotalCharge);
311       fTransform->Slice2Sector(fCurrentSlice,fCurrentRow,thissector,thisrow);
312       fTransform->Raw2Local(xyz,thissector,thisrow,fpad,ftime);
313       if(xyz[0]==0) LOG(AliL3Log::kError,"AliL3ClustFinder","Cluster Finder")
314         <<AliL3Log::kDec<<"Zero cluster"<<ENDLOG;
315       if(fNClusters >= fMaxNClusters)
316         {
317           LOG(AliL3Log::kError,"AliL3ClustFinder::WriteClusters","Cluster Finder")
318             <<AliL3Log::kDec<<"Too many clusters"<<ENDLOG;
319           return;
320           }  
321       fSpacePointData[counter].fCharge = list[j].fTotalCharge;
322       fSpacePointData[counter].fX = xyz[0];
323       fSpacePointData[counter].fY = xyz[1];
324       fSpacePointData[counter].fZ = xyz[2];
325       fSpacePointData[counter].fPadRow = fCurrentRow;
326       fSpacePointData[counter].fXYErr = fXYErr;
327       fSpacePointData[counter].fZErr = fZErr;
328       fSpacePointData[counter].fID = counter
329         +((fCurrentSlice&0x7f)<<25)+((fCurrentPatch&0x7)<<22);//uli
330       
331       fNClusters++;
332       counter++;
333
334     }
335
336 }