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