35c6950aac18a48d027331a6e286763614c674bf
[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   fcalcerr = kTRUE;
33 }
34
35 AliL3ClustFinderNew::~AliL3ClustFinderNew()
36 {
37 }
38
39 void AliL3ClustFinderNew::InitSlice(Int_t slice,Int_t patch,Int_t firstrow, Int_t lastrow,Int_t nmaxpoints)
40 {
41   fNClusters = 0;
42   fMaxNClusters = nmaxpoints;
43   fCurrentSlice = slice;
44   fCurrentPatch = patch;
45   fFirstRow = firstrow;
46   fLastRow = lastrow;
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   fSpacePointData = pt;
60 }
61
62 void AliL3ClustFinderNew::Read(UInt_t ndigits,AliL3DigitRowData *ptr)
63 {
64   fNDigitRowData = ndigits;
65   fDigitRowData = ptr;
66 }
67
68 void AliL3ClustFinderNew::ProcessDigits()
69 {
70   //Loop over rows, and call processrow
71   
72   
73   AliL3DigitRowData *tempPt = (AliL3DigitRowData*)fDigitRowData;
74   
75   for(Int_t i=fFirstRow; i<=fLastRow; i++)
76     {
77       fCurrentRow = i;
78       ProcessRow(tempPt);
79       Byte_t *tmp = (Byte_t*)tempPt;
80       Int_t size = sizeof(AliL3DigitRowData) + tempPt->fNDigit*sizeof(AliL3DigitData);
81       tmp += size;
82       tempPt = (AliL3DigitRowData*)tmp;
83     }
84   LOG(AliL3Log::kInformational,"AliL3ClustFinderNew::WriteClusters","Space points")
85     <<"Cluster finder found "<<fNClusters<<" clusters in slice "<<fCurrentSlice
86     <<" patch "<<fCurrentPatch<<ENDLOG;
87 }
88
89 void AliL3ClustFinderNew::ProcessRow(AliL3DigitRowData *tempPt)
90 {
91
92   UInt_t last_pad = 123456789;
93
94   ClusterData *pad1[2500]; //2 lists for internal memory=2pads
95   ClusterData *pad2[2500]; //2 lists for internal memory=2pads
96   ClusterData clusterlist[5000]; //Clusterlist
97
98   ClusterData **currentPt;  //List of pointers to the current pad
99   ClusterData **previousPt; //List of pointers to the previous pad
100   currentPt = pad2;
101   previousPt = pad1;
102   UInt_t n_previous=0,n_current=0,n_total=0;
103
104   //Loop over sequences of this row:
105   for(UInt_t bin=0; bin<tempPt->fNDigit; bin++)
106     {
107       
108       AliL3DigitData *data = tempPt->fDigitData;
109       if(data[bin].fPad != last_pad)
110         {
111           //This is a new pad
112           
113           //Switch:
114           if(currentPt == pad2)
115             {
116               currentPt = pad1;
117               previousPt = pad2;
118             }
119           else 
120             {
121               currentPt = pad2;
122               previousPt = pad1;
123             }
124           n_previous = n_current;
125           n_current = 0;
126           if(bin[data].fPad != last_pad+1)
127             {
128               //this happens if there is a pad with no signal.
129               n_previous = n_current = 0;
130             }
131           last_pad = data[bin].fPad;
132         }
133
134       Bool_t new_cluster = kTRUE;
135       UInt_t seq_charge=0,seq_average=0,seq_error=0;
136       UInt_t last_charge=0,last_was_falling=0;
137       Int_t new_bin=-1;
138
139       if(fDeconvTime)
140         {
141         redo: //This is a goto.
142           if(new_bin > -1)
143             {
144               bin = new_bin;
145               new_bin = -1;
146             }
147           
148           last_charge=0;
149           last_was_falling = 0;
150         }
151       
152       while(1) //Loop over current sequence
153         {
154           if(data[bin].fTime >= AliL3Transform::GetNTimeBins())
155             {
156               LOG(AliL3Log::kFatal,"AliL3ClustFinderNew::ProcessRow","Digits")
157                 <<"Timebin out of range "<<(Int_t)data[bin].fTime<<ENDLOG;
158               break;
159             }
160
161           //Get the current ADC-value
162           UInt_t charge = data[bin].fCharge;
163           
164           if(fDeconvTime)
165             {
166               //Check if the last pixel in the sequence is smaller than this
167               if(charge > last_charge)
168                 {
169                   if(last_was_falling)
170                     {
171                       new_bin = bin;
172                       break;
173                     }
174                 }
175               else last_was_falling = 1; //last pixel was larger than this
176               last_charge = charge;
177             }
178           
179           //Sum the total charge of this sequence
180           seq_charge += charge;
181           seq_average += data[bin].fTime*charge;
182           seq_error += data[bin].fTime*data[bin].fTime*charge;
183
184           //Check where to stop:
185           if(data[bin+1].fPad != data[bin].fPad) //new pad
186             break; 
187           if(data[bin+1].fTime != data[bin].fTime+1) //end of sequence
188             break;
189           
190           bin++;
191         }//end loop over sequence
192       
193       //Calculate mean of sequence:
194       Int_t seq_mean=0;
195       if(seq_charge)
196         seq_mean = seq_average/seq_charge;
197       else
198         {
199           LOG(AliL3Log::kFatal,"AliL3ClustFinderNew::ProcessRow","Data")
200             <<"Error in data given to the cluster finder"<<ENDLOG;
201           seq_mean = 1;
202           seq_charge = 1;
203         }
204       
205       //Calculate mean in pad direction:
206       Int_t pad_mean = seq_charge*data[bin].fPad;
207       Int_t pad_error = data[bin].fPad*pad_mean;
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->fPad2 += pad_error;
243               local->fTime += seq_average;
244               local->fTime2 += seq_error;
245               local->fMean = seq_mean;
246               local->fFlags++; //means we have more than one pad 
247               
248               currentPt[n_current] = local;
249               n_current++;
250               
251               break;
252             } //Checking for match at previous pad
253         } //Loop over results on previous pad.
254       
255       if(new_cluster)
256         {
257           //Start a new cluster. Add it to the clusterlist, and update
258           //the list of pointers to clusters in current pad.
259           //current pad will be previous pad on next pad.
260
261           //Add to the clusterlist:
262           ClusterData *tmp = &clusterlist[n_total];
263           tmp->fTotalCharge = seq_charge;
264           tmp->fPad = pad_mean;
265           tmp->fPad2 = pad_error;
266           tmp->fTime = seq_average;
267           tmp->fTime2 = seq_error;
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 fpad2=fXYErr;
302       Float_t ftime =(Float_t)list[j].fTime /(Float_t)list[j].fTotalCharge;
303       Float_t ftime2=fZErr;
304
305       if(fcalcerr) {
306         fpad2=(Float_t)list[j].fPad2/(Float_t)list[j].fTotalCharge - fpad*fpad;
307         fpad2 = sqrt(fpad2);
308         ftime2=(Float_t)list[j].fTime2/(Float_t)list[j].fTotalCharge - ftime*ftime;
309         ftime2 = sqrt(ftime2); 
310       }
311        
312       if(fstdout==kTRUE)
313         cout<<"WriteCluster: padrow "<<fCurrentRow<<" pad "<<fpad << "+-"<<fpad2<<" time "<<ftime<<"+-"<<ftime2<<" charge "<<list[j].fTotalCharge<<endl;
314
315       AliL3Transform::Slice2Sector(fCurrentSlice,fCurrentRow,thissector,thisrow);
316       AliL3Transform::Raw2Local(xyz,thissector,thisrow,fpad,ftime);
317       if(xyz[0]==0) LOG(AliL3Log::kError,"AliL3ClustFinder","Cluster Finder")
318         <<AliL3Log::kDec<<"Zero cluster"<<ENDLOG;
319       if(fNClusters >= fMaxNClusters)
320         {
321           LOG(AliL3Log::kError,"AliL3ClustFinder::WriteClusters","Cluster Finder")
322             <<AliL3Log::kDec<<"Too many clusters"<<ENDLOG;
323           return;
324         }  
325       fSpacePointData[counter].fCharge = list[j].fTotalCharge;
326       fSpacePointData[counter].fX = xyz[0];
327       fSpacePointData[counter].fY = xyz[1];
328       fSpacePointData[counter].fZ = xyz[2];
329       fSpacePointData[counter].fPadRow = fCurrentRow;
330       fSpacePointData[counter].fXYErr = fpad2;
331       fSpacePointData[counter].fZErr  = ftime2;
332       fSpacePointData[counter].fID = counter
333         +((fCurrentSlice&0x7f)<<25)+((fCurrentPatch&0x7)<<22);//Uli
334 #ifdef do_mc
335       Int_t trackID[3];
336       GetTrackID((Int_t)rint(fpad),(Int_t)rint(ftime),trackID);
337       //cout<<"padrow "<<fCurrentRow<<" pad "<<(Int_t)rint(fpad)<<" time "<<(Int_t)rint(ftime)<<" Trackid "<<trackID[0]<<endl;
338       fSpacePointData[counter].fTrackID[0] = trackID[0];
339       fSpacePointData[counter].fTrackID[1] = trackID[1];
340       fSpacePointData[counter].fTrackID[2] = trackID[2];
341 #endif
342       fNClusters++;
343       counter++;
344     }
345 }
346
347 #ifdef do_mc
348 void AliL3ClustFinderNew::GetTrackID(Int_t pad,Int_t time,Int_t *trackID)
349 {
350   AliL3DigitRowData *rowPt = (AliL3DigitRowData*)fDigitRowData;
351   
352   trackID[0]=trackID[1]=trackID[2]=-2;
353   //cout<<"Looking for pad "<<pad<<" time "<<time<<endl;
354   for(Int_t i=fFirstRow; i<=fLastRow; i++)
355     {
356       if(rowPt->fRow < (UInt_t)fCurrentRow)
357         {
358           AliL3MemHandler::UpdateRowPointer(rowPt);
359           continue;
360         }
361       AliL3DigitData *digPt = (AliL3DigitData*)rowPt->fDigitData;
362       for(UInt_t j=0; j<rowPt->fNDigit; j++)
363         {
364           Int_t cpad = digPt[j].fPad;
365           Int_t ctime = digPt[j].fTime;
366           if(cpad != pad) continue;
367           if(ctime != time) continue;
368           //if(cpad != pad && ctime != ctime) continue;
369           //cout<<"Reading row "<<fCurrentRow<<" pad "<<cpad<<" time "<<ctime<<" trackID "<<digPt[j].fTrackID[0]<<endl;
370           trackID[0] = digPt[j].fTrackID[0];
371           trackID[1] = digPt[j].fTrackID[1];
372           trackID[2] = digPt[j].fTrackID[2];
373           break;
374           //cout<<"Reading trackID "<<trackID[0]<<endl;
375         }
376       break;
377     }
378   
379 }
380 #endif