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