Solved bug to not merge clusters on the same pad.
[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
b62fa469 15/** \class AliL3ClustFinderNew
16//<pre>
b661165c 17//_____________________________________________________________
18// AliL3ClustFinderNew
19//
20// The current cluster finder for HLT
21// Based on STAR L3
b62fa469 22//</pre> */
b661165c 23
46fbeb8a 24ClassImp(AliL3ClustFinderNew)
25
26AliL3ClustFinderNew::AliL3ClustFinderNew()
27{
28 fMatch = 4;
77f08cd3 29 fThreshold = 10;
30 fXYErr = 0.2;
31 fZErr = 0.3;
99a6d5c9 32 fDeconvPad = kTRUE;
2bbf57d8 33 fDeconvTime = kTRUE;
77f08cd3 34 fstdout = kFALSE;
c3dd27a3 35 fcalcerr = kTRUE;
46fbeb8a 36}
37
46fbeb8a 38AliL3ClustFinderNew::~AliL3ClustFinderNew()
39{
46fbeb8a 40}
41
42void 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
52void 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
60void AliL3ClustFinderNew::SetOutputArray(AliL3SpacePointData *pt)
61{
46fbeb8a 62 fSpacePointData = pt;
63}
64
46fbeb8a 65void AliL3ClustFinderNew::Read(UInt_t ndigits,AliL3DigitRowData *ptr)
66{
67 fNDigitRowData = ndigits;
68 fDigitRowData = ptr;
46fbeb8a 69}
70
71void 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 }
355debe1 87 LOG(AliL3Log::kInformational,"AliL3ClustFinderNew::WriteClusters","Space points")
88 <<"Cluster finder found "<<fNClusters<<" clusters in slice "<<fCurrentSlice
89 <<" patch "<<fCurrentPatch<<ENDLOG;
46fbeb8a 90}
91
46fbeb8a 92void 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
77f08cd3 101 ClusterData **currentPt; //List of pointers to the current pad
46fbeb8a 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 {
46fbeb8a 110 AliL3DigitData *data = tempPt->fDigitData;
111 if(data[bin].fPad != last_pad)
112 {
113 //This is a new pad
46fbeb8a 114
115 //Switch:
116 if(currentPt == pad2)
117 {
118 currentPt = pad1;
119 previousPt = pad2;
46fbeb8a 120 }
121 else
122 {
123 currentPt = pad2;
124 previousPt = pad1;
125 }
126 n_previous = n_current;
127 n_current = 0;
99a6d5c9 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;
46fbeb8a 134 }
135
136 Bool_t new_cluster = kTRUE;
c3dd27a3 137 UInt_t seq_charge=0,seq_average=0,seq_error=0;
99a6d5c9 138 UInt_t last_charge=0,last_was_falling=0;
139 Int_t new_bin=-1;
77f08cd3 140
99a6d5c9 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
46fbeb8a 154 while(1) //Loop over current sequence
155 {
494fad94 156 if(data[bin].fTime >= AliL3Transform::GetNTimeBins())
adb0e845 157 {
158 LOG(AliL3Log::kFatal,"AliL3ClustFinderNew::ProcessRow","Digits")
159 <<"Timebin out of range "<<(Int_t)data[bin].fTime<<ENDLOG;
160 break;
161 }
162
46fbeb8a 163 //Get the current ADC-value
164 UInt_t charge = data[bin].fCharge;
99a6d5c9 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
46fbeb8a 181 //Sum the total charge of this sequence
182 seq_charge += charge;
46fbeb8a 183 seq_average += data[bin].fTime*charge;
c3dd27a3 184 seq_error += data[bin].fTime*data[bin].fTime*charge;
185
46fbeb8a 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
99a6d5c9 200 {
adb0e845 201 LOG(AliL3Log::kFatal,"AliL3ClustFinderNew::ProcessRow","Data")
202 <<"Error in data given to the cluster finder"<<ENDLOG;
99a6d5c9 203 seq_mean = 1;
204 seq_charge = 1;
205 }
46fbeb8a 206
207 //Calculate mean in pad direction:
208 Int_t pad_mean = seq_charge*data[bin].fPad;
c3dd27a3 209 Int_t pad_error = data[bin].fPad*pad_mean;
b62fa469 210
46fbeb8a 211 //Compare with results on previous pad:
212 for(UInt_t p=0; p<n_previous; p++)
213 {
b62fa469 214 //dont merge sequences on the same pad twice
215 if(previousPt[p]->fLastMergedPad==data[bin].fPad) continue;
216
46fbeb8a 217 Int_t difference = seq_mean - previousPt[p]->fMean;
b62fa469 218 if(difference < -fMatch) break;
46fbeb8a 219
c3dd27a3 220 if(difference <= fMatch) //There is a match here!!
46fbeb8a 221 {
99a6d5c9 222 ClusterData *local = previousPt[p];
b62fa469 223
99a6d5c9 224 if(fDeconvPad)
225 {
226 if(seq_charge > local->fLastCharge)
227 {
c3dd27a3 228 if(local->fChargeFalling) //The previous pad was falling
99a6d5c9 229 {
c3dd27a3 230 break; //create a new cluster
99a6d5c9 231 }
232 }
233 else
234 local->fChargeFalling = 1;
235 local->fLastCharge = seq_charge;
236 }
237
46fbeb8a 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:
99a6d5c9 242 local->fTotalCharge += seq_charge;
243 local->fPad += pad_mean;
c3dd27a3 244 local->fPad2 += pad_error;
99a6d5c9 245 local->fTime += seq_average;
c3dd27a3 246 local->fTime2 += seq_error;
99a6d5c9 247 local->fMean = seq_mean;
248 local->fFlags++; //means we have more than one pad
b62fa469 249 local->fLastMergedPad = data[bin].fPad;
250
99a6d5c9 251 currentPt[n_current] = local;
46fbeb8a 252 n_current++;
253
46fbeb8a 254 break;
c3dd27a3 255 } //Checking for match at previous pad
256 } //Loop over results on previous pad.
46fbeb8a 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.
99a6d5c9 262 //current pad will be previous pad on next pad.
46fbeb8a 263
264 //Add to the clusterlist:
265 ClusterData *tmp = &clusterlist[n_total];
266 tmp->fTotalCharge = seq_charge;
267 tmp->fPad = pad_mean;
c3dd27a3 268 tmp->fPad2 = pad_error;
46fbeb8a 269 tmp->fTime = seq_average;
c3dd27a3 270 tmp->fTime2 = seq_error;
46fbeb8a 271 tmp->fMean = seq_mean;
272 tmp->fFlags = 0; //flags for 1 pad clusters
b62fa469 273 tmp->fLastMergedPad = data[bin].fPad;
99a6d5c9 274 if(fDeconvPad)
275 {
276 tmp->fChargeFalling = 0;
277 tmp->fLastCharge = seq_charge;
278 }
279
46fbeb8a 280 //Update list of pointers to previous pad:
281 currentPt[n_current] = &clusterlist[n_total];
282 n_total++;
283 n_current++;
284 }
b62fa469 285
99a6d5c9 286 if(fDeconvTime)
287 if(new_bin >= 0) goto redo;
46fbeb8a 288 }//Loop over digits on this padrow
289
290 WriteClusters(n_total,clusterlist);
46fbeb8a 291}
292
293void 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
d3b44290 302
303 Float_t xyz[3];
c3dd27a3 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;
eeddc64d 308
c3dd27a3 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
77f08cd3 316 if(fstdout==kTRUE)
b62fa469 317 cout<<"WriteCluster: padrow "<<fCurrentRow<<" pad "<<fpad << " +- "<<fpad2<<" time "<<ftime<<" +- "<<ftime2<<" charge "<<list[j].fTotalCharge<<endl;
d3b44290 318
494fad94 319 AliL3Transform::Slice2Sector(fCurrentSlice,fCurrentRow,thissector,thisrow);
320 AliL3Transform::Raw2Local(xyz,thissector,thisrow,fpad,ftime);
46fbeb8a 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;
d3b44290 328 }
95a00d93 329 fSpacePointData[counter].fCharge = list[j].fTotalCharge;
46fbeb8a 330 fSpacePointData[counter].fX = xyz[0];
331 fSpacePointData[counter].fY = xyz[1];
332 fSpacePointData[counter].fZ = xyz[2];
333 fSpacePointData[counter].fPadRow = fCurrentRow;
c3dd27a3 334 fSpacePointData[counter].fXYErr = fpad2;
335 fSpacePointData[counter].fZErr = ftime2;
46fbeb8a 336 fSpacePointData[counter].fID = counter
d3b44290 337 +((fCurrentSlice&0x7f)<<25)+((fCurrentPatch&0x7)<<22);//Uli
10f815d9 338#ifdef do_mc
339 Int_t trackID[3];
340 GetTrackID((Int_t)rint(fpad),(Int_t)rint(ftime),trackID);
6ab8a12a 341 //cout<<"padrow "<<fCurrentRow<<" pad "<<(Int_t)rint(fpad)<<" time "<<(Int_t)rint(ftime)<<" Trackid "<<trackID[0]<<endl;
10f815d9 342 fSpacePointData[counter].fTrackID[0] = trackID[0];
343 fSpacePointData[counter].fTrackID[1] = trackID[1];
344 fSpacePointData[counter].fTrackID[2] = trackID[2];
345#endif
46fbeb8a 346 fNClusters++;
347 counter++;
46fbeb8a 348 }
46fbeb8a 349}
10f815d9 350
351#ifdef do_mc
352void AliL3ClustFinderNew::GetTrackID(Int_t pad,Int_t time,Int_t *trackID)
353{
354 AliL3DigitRowData *rowPt = (AliL3DigitRowData*)fDigitRowData;
355
355debe1 356 trackID[0]=trackID[1]=trackID[2]=-2;
10f815d9 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;
355debe1 370 if(cpad != pad) continue;
371 if(ctime != time) continue;
372 //if(cpad != pad && ctime != ctime) continue;
10f815d9 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];
355debe1 377 break;
10f815d9 378 //cout<<"Reading trackID "<<trackID[0]<<endl;
379 }
380 break;
381 }
382
383}
384#endif