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