]> git.uio.no Git - u/mrichter/AliRoot.git/blame - HLT/src/AliL3ClustFinderNew.cxx
Added the possibility to save the particle id's through the chain, if detailed effici...
[u/mrichter/AliRoot.git] / HLT / src / AliL3ClustFinderNew.cxx
CommitLineData
4ec96209 1//$Id$
2
b661165c 3// Author: Anders Vestbo <mailto:vestbo@fi.uib.no>
4//*-- Copyright &copy ASV
46fbeb8a 5
6#include <iostream.h>
10f815d9 7#include <math.h>
46fbeb8a 8#include "AliL3Logging.h"
9#include "AliL3ClustFinderNew.h"
10#include "AliL3DigitData.h"
11#include "AliL3Transform.h"
12#include "AliL3SpacePointData.h"
10f815d9 13#include "AliL3MemHandler.h"
46fbeb8a 14
b661165c 15//_____________________________________________________________
16// AliL3ClustFinderNew
17//
18// The current cluster finder for HLT
19// Based on STAR L3
20
46fbeb8a 21ClassImp(AliL3ClustFinderNew)
22
23AliL3ClustFinderNew::AliL3ClustFinderNew()
24{
25 fMatch = 4;
26 fThreshold =10;
99a6d5c9 27 fDeconvPad = kTRUE;
2bbf57d8 28 fDeconvTime = kTRUE;
46fbeb8a 29}
30
31AliL3ClustFinderNew::AliL3ClustFinderNew(AliL3Transform *transform)
32{
33 fTransform = transform;
34 fMatch = 4;
35 fThreshold =10;
99a6d5c9 36 fDeconvTime = kTRUE;
37 fDeconvPad = kTRUE;
46fbeb8a 38}
39
40AliL3ClustFinderNew::~AliL3ClustFinderNew()
41{
42
43
44}
45
46void 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;
99a6d5c9 54 fDeconvTime = kTRUE;
55 fDeconvPad = kTRUE;
46fbeb8a 56}
57
58void 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
66void AliL3ClustFinderNew::SetOutputArray(AliL3SpacePointData *pt)
67{
68
69 fSpacePointData = pt;
70}
71
72
73void AliL3ClustFinderNew::Read(UInt_t ndigits,AliL3DigitRowData *ptr)
74{
75 fNDigitRowData = ndigits;
76 fDigitRowData = ptr;
77
78}
79
80void 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
100void 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
46fbeb8a 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;
99a6d5c9 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;
46fbeb8a 144 }
145
146 Bool_t new_cluster = kTRUE;
147
148 UInt_t seq_charge=0,seq_average=0;
149
99a6d5c9 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
46fbeb8a 166 while(1) //Loop over current sequence
167 {
adb0e845 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
46fbeb8a 175 //Get the current ADC-value
176 UInt_t charge = data[bin].fCharge;
99a6d5c9 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
46fbeb8a 193 //Sum the total charge of this sequence
194 seq_charge += charge;
46fbeb8a 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
99a6d5c9 211 {
adb0e845 212 LOG(AliL3Log::kFatal,"AliL3ClustFinderNew::ProcessRow","Data")
213 <<"Error in data given to the cluster finder"<<ENDLOG;
99a6d5c9 214 seq_mean = 1;
215 seq_charge = 1;
216 }
46fbeb8a 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 {
99a6d5c9 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
46fbeb8a 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:
99a6d5c9 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;
46fbeb8a 259 n_current++;
260
46fbeb8a 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.
99a6d5c9 269 //current pad will be previous pad on next pad.
46fbeb8a 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
99a6d5c9 278 if(fDeconvPad)
279 {
280 tmp->fChargeFalling = 0;
281 tmp->fLastCharge = seq_charge;
282 }
283
46fbeb8a 284 //Update list of pointers to previous pad:
285 currentPt[n_current] = &clusterlist[n_total];
286 n_total++;
287 n_current++;
288 }
99a6d5c9 289 if(fDeconvTime)
290 if(new_bin >= 0) goto redo;
291
46fbeb8a 292 }//Loop over digits on this padrow
293
294 WriteClusters(n_total,clusterlist);
295
296
297}
298
299void 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 }
95a00d93 323 fSpacePointData[counter].fCharge = list[j].fTotalCharge;
46fbeb8a 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
10f815d9 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
46fbeb8a 339 fNClusters++;
340 counter++;
341
342 }
343
344}
10f815d9 345
346#ifdef do_mc
347void 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