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