Default match should be 2
[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
118c26c3 6
7#include "AliL3StandardIncludes.h"
8
46fbeb8a 9#include "AliL3Logging.h"
10#include "AliL3ClustFinderNew.h"
11#include "AliL3DigitData.h"
12#include "AliL3Transform.h"
13#include "AliL3SpacePointData.h"
10f815d9 14#include "AliL3MemHandler.h"
46fbeb8a 15
118c26c3 16#if GCCVERSION == 3
17using namespace std;
18#endif
19
b62fa469 20/** \class AliL3ClustFinderNew
21//<pre>
b661165c 22//_____________________________________________________________
23// AliL3ClustFinderNew
24//
25// The current cluster finder for HLT
26// Based on STAR L3
b62fa469 27//</pre> */
b661165c 28
46fbeb8a 29ClassImp(AliL3ClustFinderNew)
30
31AliL3ClustFinderNew::AliL3ClustFinderNew()
32{
9953e05f 33 fMatch = 2;
77f08cd3 34 fThreshold = 10;
35 fXYErr = 0.2;
36 fZErr = 0.3;
99a6d5c9 37 fDeconvPad = kTRUE;
2bbf57d8 38 fDeconvTime = kTRUE;
77f08cd3 39 fstdout = kFALSE;
c3dd27a3 40 fcalcerr = kTRUE;
46fbeb8a 41}
42
46fbeb8a 43AliL3ClustFinderNew::~AliL3ClustFinderNew()
44{
46fbeb8a 45}
46
47void 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
57void 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
65void AliL3ClustFinderNew::SetOutputArray(AliL3SpacePointData *pt)
66{
46fbeb8a 67 fSpacePointData = pt;
68}
69
46fbeb8a 70void AliL3ClustFinderNew::Read(UInt_t ndigits,AliL3DigitRowData *ptr)
71{
72 fNDigitRowData = ndigits;
73 fDigitRowData = ptr;
46fbeb8a 74}
75
76void 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 }
355debe1 92 LOG(AliL3Log::kInformational,"AliL3ClustFinderNew::WriteClusters","Space points")
93 <<"Cluster finder found "<<fNClusters<<" clusters in slice "<<fCurrentSlice
94 <<" patch "<<fCurrentPatch<<ENDLOG;
46fbeb8a 95}
96
46fbeb8a 97void 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
77f08cd3 106 ClusterData **currentPt; //List of pointers to the current pad
46fbeb8a 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 {
46fbeb8a 115 AliL3DigitData *data = tempPt->fDigitData;
116 if(data[bin].fPad != last_pad)
117 {
118 //This is a new pad
46fbeb8a 119
120 //Switch:
121 if(currentPt == pad2)
122 {
123 currentPt = pad1;
124 previousPt = pad2;
46fbeb8a 125 }
126 else
127 {
128 currentPt = pad2;
129 previousPt = pad1;
130 }
131 n_previous = n_current;
132 n_current = 0;
99a6d5c9 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;
46fbeb8a 139 }
140
141 Bool_t new_cluster = kTRUE;
c3dd27a3 142 UInt_t seq_charge=0,seq_average=0,seq_error=0;
99a6d5c9 143 UInt_t last_charge=0,last_was_falling=0;
144 Int_t new_bin=-1;
77f08cd3 145
99a6d5c9 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
46fbeb8a 159 while(1) //Loop over current sequence
160 {
494fad94 161 if(data[bin].fTime >= AliL3Transform::GetNTimeBins())
adb0e845 162 {
163 LOG(AliL3Log::kFatal,"AliL3ClustFinderNew::ProcessRow","Digits")
164 <<"Timebin out of range "<<(Int_t)data[bin].fTime<<ENDLOG;
165 break;
166 }
167
46fbeb8a 168 //Get the current ADC-value
169 UInt_t charge = data[bin].fCharge;
99a6d5c9 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
46fbeb8a 186 //Sum the total charge of this sequence
187 seq_charge += charge;
46fbeb8a 188 seq_average += data[bin].fTime*charge;
c3dd27a3 189 seq_error += data[bin].fTime*data[bin].fTime*charge;
190
46fbeb8a 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
99a6d5c9 205 {
adb0e845 206 LOG(AliL3Log::kFatal,"AliL3ClustFinderNew::ProcessRow","Data")
207 <<"Error in data given to the cluster finder"<<ENDLOG;
99a6d5c9 208 seq_mean = 1;
209 seq_charge = 1;
210 }
46fbeb8a 211
212 //Calculate mean in pad direction:
213 Int_t pad_mean = seq_charge*data[bin].fPad;
c3dd27a3 214 Int_t pad_error = data[bin].fPad*pad_mean;
b62fa469 215
46fbeb8a 216 //Compare with results on previous pad:
217 for(UInt_t p=0; p<n_previous; p++)
218 {
b62fa469 219 //dont merge sequences on the same pad twice
220 if(previousPt[p]->fLastMergedPad==data[bin].fPad) continue;
221
46fbeb8a 222 Int_t difference = seq_mean - previousPt[p]->fMean;
b62fa469 223 if(difference < -fMatch) break;
46fbeb8a 224
c3dd27a3 225 if(difference <= fMatch) //There is a match here!!
46fbeb8a 226 {
99a6d5c9 227 ClusterData *local = previousPt[p];
b62fa469 228
99a6d5c9 229 if(fDeconvPad)
230 {
231 if(seq_charge > local->fLastCharge)
232 {
c3dd27a3 233 if(local->fChargeFalling) //The previous pad was falling
99a6d5c9 234 {
c3dd27a3 235 break; //create a new cluster
99a6d5c9 236 }
237 }
238 else
239 local->fChargeFalling = 1;
240 local->fLastCharge = seq_charge;
241 }
242
46fbeb8a 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:
99a6d5c9 247 local->fTotalCharge += seq_charge;
248 local->fPad += pad_mean;
c3dd27a3 249 local->fPad2 += pad_error;
99a6d5c9 250 local->fTime += seq_average;
c3dd27a3 251 local->fTime2 += seq_error;
99a6d5c9 252 local->fMean = seq_mean;
253 local->fFlags++; //means we have more than one pad
b62fa469 254 local->fLastMergedPad = data[bin].fPad;
255
99a6d5c9 256 currentPt[n_current] = local;
46fbeb8a 257 n_current++;
258
46fbeb8a 259 break;
c3dd27a3 260 } //Checking for match at previous pad
261 } //Loop over results on previous pad.
46fbeb8a 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.
99a6d5c9 267 //current pad will be previous pad on next pad.
46fbeb8a 268
269 //Add to the clusterlist:
270 ClusterData *tmp = &clusterlist[n_total];
271 tmp->fTotalCharge = seq_charge;
272 tmp->fPad = pad_mean;
c3dd27a3 273 tmp->fPad2 = pad_error;
46fbeb8a 274 tmp->fTime = seq_average;
c3dd27a3 275 tmp->fTime2 = seq_error;
46fbeb8a 276 tmp->fMean = seq_mean;
277 tmp->fFlags = 0; //flags for 1 pad clusters
b62fa469 278 tmp->fLastMergedPad = data[bin].fPad;
99a6d5c9 279 if(fDeconvPad)
280 {
281 tmp->fChargeFalling = 0;
282 tmp->fLastCharge = seq_charge;
283 }
284
46fbeb8a 285 //Update list of pointers to previous pad:
286 currentPt[n_current] = &clusterlist[n_total];
287 n_total++;
288 n_current++;
289 }
b62fa469 290
99a6d5c9 291 if(fDeconvTime)
292 if(new_bin >= 0) goto redo;
46fbeb8a 293 }//Loop over digits on this padrow
294
295 WriteClusters(n_total,clusterlist);
46fbeb8a 296}
297
298void 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
d3b44290 307
308 Float_t xyz[3];
c3dd27a3 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;
eeddc64d 313
c3dd27a3 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
77f08cd3 321 if(fstdout==kTRUE)
b62fa469 322 cout<<"WriteCluster: padrow "<<fCurrentRow<<" pad "<<fpad << " +- "<<fpad2<<" time "<<ftime<<" +- "<<ftime2<<" charge "<<list[j].fTotalCharge<<endl;
d3b44290 323
494fad94 324 AliL3Transform::Slice2Sector(fCurrentSlice,fCurrentRow,thissector,thisrow);
325 AliL3Transform::Raw2Local(xyz,thissector,thisrow,fpad,ftime);
46fbeb8a 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;
d3b44290 333 }
95a00d93 334 fSpacePointData[counter].fCharge = list[j].fTotalCharge;
46fbeb8a 335 fSpacePointData[counter].fX = xyz[0];
336 fSpacePointData[counter].fY = xyz[1];
337 fSpacePointData[counter].fZ = xyz[2];
338 fSpacePointData[counter].fPadRow = fCurrentRow;
c3dd27a3 339 fSpacePointData[counter].fXYErr = fpad2;
340 fSpacePointData[counter].fZErr = ftime2;
46fbeb8a 341 fSpacePointData[counter].fID = counter
d3b44290 342 +((fCurrentSlice&0x7f)<<25)+((fCurrentPatch&0x7)<<22);//Uli
10f815d9 343#ifdef do_mc
344 Int_t trackID[3];
345 GetTrackID((Int_t)rint(fpad),(Int_t)rint(ftime),trackID);
6ab8a12a 346 //cout<<"padrow "<<fCurrentRow<<" pad "<<(Int_t)rint(fpad)<<" time "<<(Int_t)rint(ftime)<<" Trackid "<<trackID[0]<<endl;
10f815d9 347 fSpacePointData[counter].fTrackID[0] = trackID[0];
348 fSpacePointData[counter].fTrackID[1] = trackID[1];
349 fSpacePointData[counter].fTrackID[2] = trackID[2];
350#endif
9953e05f 351
46fbeb8a 352 fNClusters++;
353 counter++;
46fbeb8a 354 }
46fbeb8a 355}
10f815d9 356
357#ifdef do_mc
358void AliL3ClustFinderNew::GetTrackID(Int_t pad,Int_t time,Int_t *trackID)
359{
360 AliL3DigitRowData *rowPt = (AliL3DigitRowData*)fDigitRowData;
361
355debe1 362 trackID[0]=trackID[1]=trackID[2]=-2;
10f815d9 363 //cout<<"Looking for pad "<<pad<<" time "<<time<<endl;
364 for(Int_t i=fFirstRow; i<=fLastRow; i++)
365 {
366 if(rowPt->fRow < (UInt_t)fCurrentRow)
367 {
368 AliL3MemHandler::UpdateRowPointer(rowPt);
369 continue;
370 }
371 AliL3DigitData *digPt = (AliL3DigitData*)rowPt->fDigitData;
372 for(UInt_t j=0; j<rowPt->fNDigit; j++)
373 {
374 Int_t cpad = digPt[j].fPad;
375 Int_t ctime = digPt[j].fTime;
355debe1 376 if(cpad != pad) continue;
377 if(ctime != time) continue;
378 //if(cpad != pad && ctime != ctime) continue;
10f815d9 379 //cout<<"Reading row "<<fCurrentRow<<" pad "<<cpad<<" time "<<ctime<<" trackID "<<digPt[j].fTrackID[0]<<endl;
380 trackID[0] = digPt[j].fTrackID[0];
381 trackID[1] = digPt[j].fTrackID[1];
382 trackID[2] = digPt[j].fTrackID[2];
355debe1 383 break;
10f815d9 384 //cout<<"Reading trackID "<<trackID[0]<<endl;
385 }
386 break;
387 }
388
389}
390#endif