Solved bug to not merge clusters on the same pad.
[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 /** \class AliL3ClustFinderNew
16 //<pre>
17 //_____________________________________________________________
18 // AliL3ClustFinderNew
19 //
20 // The current cluster finder for HLT
21 // Based on STAR L3
22 //</pre> */
23
24 ClassImp(AliL3ClustFinderNew)
25
26 AliL3ClustFinderNew::AliL3ClustFinderNew()
27 {
28   fMatch = 4;
29   fThreshold = 10;
30   fXYErr = 0.2;
31   fZErr = 0.3;
32   fDeconvPad = kTRUE;
33   fDeconvTime = kTRUE;
34   fstdout = kFALSE;
35   fcalcerr = kTRUE;
36 }
37
38 AliL3ClustFinderNew::~AliL3ClustFinderNew()
39 {
40 }
41
42 void AliL3ClustFinderNew::InitSlice(Int_t slice,Int_t patch,Int_t firstrow, Int_t lastrow,Int_t nmaxpoints)
43 {
44   fNClusters = 0;
45   fMaxNClusters = nmaxpoints;
46   fCurrentSlice = slice;
47   fCurrentPatch = patch;
48   fFirstRow = firstrow;
49   fLastRow = lastrow;
50 }
51
52 void AliL3ClustFinderNew::InitSlice(Int_t slice,Int_t patch,Int_t nmaxpoints)
53 {
54   fNClusters = 0;
55   fMaxNClusters = nmaxpoints;
56   fCurrentSlice = slice;
57   fCurrentPatch = patch;
58 }
59
60 void AliL3ClustFinderNew::SetOutputArray(AliL3SpacePointData *pt)
61 {
62   fSpacePointData = pt;
63 }
64
65 void AliL3ClustFinderNew::Read(UInt_t ndigits,AliL3DigitRowData *ptr)
66 {
67   fNDigitRowData = ndigits;
68   fDigitRowData = ptr;
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       AliL3DigitData *data = tempPt->fDigitData;
111       if(data[bin].fPad != last_pad)
112         {
113           //This is a new pad
114           
115           //Switch:
116           if(currentPt == pad2)
117             {
118               currentPt = pad1;
119               previousPt = pad2;
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,seq_error=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           seq_error += data[bin].fTime*data[bin].fTime*charge;
185
186           //Check where to stop:
187           if(data[bin+1].fPad != data[bin].fPad) //new pad
188             break; 
189           if(data[bin+1].fTime != data[bin].fTime+1) //end of sequence
190             break;
191           
192           bin++;
193         }//end loop over sequence
194       
195       //Calculate mean of sequence:
196       Int_t seq_mean=0;
197       if(seq_charge)
198         seq_mean = seq_average/seq_charge;
199       else
200         {
201           LOG(AliL3Log::kFatal,"AliL3ClustFinderNew::ProcessRow","Data")
202             <<"Error in data given to the cluster finder"<<ENDLOG;
203           seq_mean = 1;
204           seq_charge = 1;
205         }
206       
207       //Calculate mean in pad direction:
208       Int_t pad_mean = seq_charge*data[bin].fPad;
209       Int_t pad_error = data[bin].fPad*pad_mean;
210       
211       //Compare with results on previous pad:
212       for(UInt_t p=0; p<n_previous; p++)
213         {
214           //dont merge sequences on the same pad twice
215           if(previousPt[p]->fLastMergedPad==data[bin].fPad) continue;
216
217           Int_t difference = seq_mean - previousPt[p]->fMean;
218           if(difference < -fMatch) break;
219
220           if(difference <= fMatch) //There is a match here!!
221             {
222               ClusterData *local = previousPt[p];
223
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               local->fTotalCharge += seq_charge;
243               local->fPad += pad_mean;
244               local->fPad2 += pad_error;
245               local->fTime += seq_average;
246               local->fTime2 += seq_error;
247               local->fMean = seq_mean;
248               local->fFlags++; //means we have more than one pad 
249               local->fLastMergedPad = data[bin].fPad;
250
251               currentPt[n_current] = local;
252               n_current++;
253               
254               break;
255             } //Checking for match at previous pad
256         } //Loop over results on previous pad.
257       
258       if(new_cluster)
259         {
260           //Start a new cluster. Add it to the clusterlist, and update
261           //the list of pointers to clusters in current pad.
262           //current pad will be previous pad on next pad.
263
264           //Add to the clusterlist:
265           ClusterData *tmp = &clusterlist[n_total];
266           tmp->fTotalCharge = seq_charge;
267           tmp->fPad = pad_mean;
268           tmp->fPad2 = pad_error;
269           tmp->fTime = seq_average;
270           tmp->fTime2 = seq_error;
271           tmp->fMean = seq_mean;
272           tmp->fFlags = 0;  //flags for 1 pad clusters
273           tmp->fLastMergedPad = data[bin].fPad;
274           if(fDeconvPad)
275             {
276               tmp->fChargeFalling = 0;
277               tmp->fLastCharge = seq_charge;
278             }
279
280           //Update list of pointers to previous pad:
281           currentPt[n_current] = &clusterlist[n_total];
282           n_total++;
283           n_current++;
284         }
285
286       if(fDeconvTime)
287         if(new_bin >= 0) goto redo;
288     }//Loop over digits on this padrow
289   
290   WriteClusters(n_total,clusterlist);
291 }
292
293 void AliL3ClustFinderNew::WriteClusters(Int_t n_clusters,ClusterData *list)
294 {
295   Int_t thisrow,thissector;
296   UInt_t counter = fNClusters;
297   
298   for(int j=0; j<n_clusters; j++)
299     {
300       if(!list[j].fFlags) continue; //discard 1 pad clusters
301       if(list[j].fTotalCharge < fThreshold) continue; //noise cluster
302
303       Float_t xyz[3];      
304       Float_t fpad =(Float_t)list[j].fPad /(Float_t)list[j].fTotalCharge;
305       Float_t fpad2=fXYErr;
306       Float_t ftime =(Float_t)list[j].fTime /(Float_t)list[j].fTotalCharge;
307       Float_t ftime2=fZErr;
308
309       if(fcalcerr) {
310         fpad2=(Float_t)list[j].fPad2/(Float_t)list[j].fTotalCharge - fpad*fpad;
311         fpad2 = sqrt(fpad2);
312         ftime2=(Float_t)list[j].fTime2/(Float_t)list[j].fTotalCharge - ftime*ftime;
313         ftime2 = sqrt(ftime2); 
314       }
315        
316       if(fstdout==kTRUE)
317         cout<<"WriteCluster: padrow "<<fCurrentRow<<" pad "<<fpad << " +- "<<fpad2<<" time "<<ftime<<" +- "<<ftime2<<" charge "<<list[j].fTotalCharge<<endl;
318
319       AliL3Transform::Slice2Sector(fCurrentSlice,fCurrentRow,thissector,thisrow);
320       AliL3Transform::Raw2Local(xyz,thissector,thisrow,fpad,ftime);
321       if(xyz[0]==0) LOG(AliL3Log::kError,"AliL3ClustFinder","Cluster Finder")
322         <<AliL3Log::kDec<<"Zero cluster"<<ENDLOG;
323       if(fNClusters >= fMaxNClusters)
324         {
325           LOG(AliL3Log::kError,"AliL3ClustFinder::WriteClusters","Cluster Finder")
326             <<AliL3Log::kDec<<"Too many clusters"<<ENDLOG;
327           return;
328         }  
329       fSpacePointData[counter].fCharge = list[j].fTotalCharge;
330       fSpacePointData[counter].fX = xyz[0];
331       fSpacePointData[counter].fY = xyz[1];
332       fSpacePointData[counter].fZ = xyz[2];
333       fSpacePointData[counter].fPadRow = fCurrentRow;
334       fSpacePointData[counter].fXYErr = fpad2;
335       fSpacePointData[counter].fZErr  = ftime2;
336       fSpacePointData[counter].fID = counter
337         +((fCurrentSlice&0x7f)<<25)+((fCurrentPatch&0x7)<<22);//Uli
338 #ifdef do_mc
339       Int_t trackID[3];
340       GetTrackID((Int_t)rint(fpad),(Int_t)rint(ftime),trackID);
341       //cout<<"padrow "<<fCurrentRow<<" pad "<<(Int_t)rint(fpad)<<" time "<<(Int_t)rint(ftime)<<" Trackid "<<trackID[0]<<endl;
342       fSpacePointData[counter].fTrackID[0] = trackID[0];
343       fSpacePointData[counter].fTrackID[1] = trackID[1];
344       fSpacePointData[counter].fTrackID[2] = trackID[2];
345 #endif
346       fNClusters++;
347       counter++;
348     }
349 }
350
351 #ifdef do_mc
352 void AliL3ClustFinderNew::GetTrackID(Int_t pad,Int_t time,Int_t *trackID)
353 {
354   AliL3DigitRowData *rowPt = (AliL3DigitRowData*)fDigitRowData;
355   
356   trackID[0]=trackID[1]=trackID[2]=-2;
357   //cout<<"Looking for pad "<<pad<<" time "<<time<<endl;
358   for(Int_t i=fFirstRow; i<=fLastRow; i++)
359     {
360       if(rowPt->fRow < (UInt_t)fCurrentRow)
361         {
362           AliL3MemHandler::UpdateRowPointer(rowPt);
363           continue;
364         }
365       AliL3DigitData *digPt = (AliL3DigitData*)rowPt->fDigitData;
366       for(UInt_t j=0; j<rowPt->fNDigit; j++)
367         {
368           Int_t cpad = digPt[j].fPad;
369           Int_t ctime = digPt[j].fTime;
370           if(cpad != pad) continue;
371           if(ctime != time) continue;
372           //if(cpad != pad && ctime != ctime) continue;
373           //cout<<"Reading row "<<fCurrentRow<<" pad "<<cpad<<" time "<<ctime<<" trackID "<<digPt[j].fTrackID[0]<<endl;
374           trackID[0] = digPt[j].fTrackID[0];
375           trackID[1] = digPt[j].fTrackID[1];
376           trackID[2] = digPt[j].fTrackID[2];
377           break;
378           //cout<<"Reading trackID "<<trackID[0]<<endl;
379         }
380       break;
381     }
382   
383 }
384 #endif