Changes done to make the Cluser Finder calculate the errors in Pad and Time direction...
[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;
77f08cd3 26 fThreshold = 10;
27 fXYErr = 0.2;
28 fZErr = 0.3;
99a6d5c9 29 fDeconvPad = kTRUE;
2bbf57d8 30 fDeconvTime = kTRUE;
77f08cd3 31 fstdout = kFALSE;
c3dd27a3 32 fcalcerr = kTRUE;
46fbeb8a 33}
34
46fbeb8a 35AliL3ClustFinderNew::~AliL3ClustFinderNew()
36{
46fbeb8a 37}
38
39void 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
49void 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
57void AliL3ClustFinderNew::SetOutputArray(AliL3SpacePointData *pt)
58{
46fbeb8a 59 fSpacePointData = pt;
60}
61
46fbeb8a 62void AliL3ClustFinderNew::Read(UInt_t ndigits,AliL3DigitRowData *ptr)
63{
64 fNDigitRowData = ndigits;
65 fDigitRowData = ptr;
46fbeb8a 66}
67
68void 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 }
355debe1 84 LOG(AliL3Log::kInformational,"AliL3ClustFinderNew::WriteClusters","Space points")
85 <<"Cluster finder found "<<fNClusters<<" clusters in slice "<<fCurrentSlice
86 <<" patch "<<fCurrentPatch<<ENDLOG;
46fbeb8a 87}
88
46fbeb8a 89void 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
77f08cd3 98 ClusterData **currentPt; //List of pointers to the current pad
46fbeb8a 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
46fbeb8a 112
113 //Switch:
114 if(currentPt == pad2)
115 {
116 currentPt = pad1;
117 previousPt = pad2;
46fbeb8a 118 }
119 else
120 {
121 currentPt = pad2;
122 previousPt = pad1;
123 }
124 n_previous = n_current;
125 n_current = 0;
99a6d5c9 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;
46fbeb8a 132 }
133
134 Bool_t new_cluster = kTRUE;
c3dd27a3 135 UInt_t seq_charge=0,seq_average=0,seq_error=0;
99a6d5c9 136 UInt_t last_charge=0,last_was_falling=0;
137 Int_t new_bin=-1;
77f08cd3 138
99a6d5c9 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
46fbeb8a 152 while(1) //Loop over current sequence
153 {
494fad94 154 if(data[bin].fTime >= AliL3Transform::GetNTimeBins())
adb0e845 155 {
156 LOG(AliL3Log::kFatal,"AliL3ClustFinderNew::ProcessRow","Digits")
157 <<"Timebin out of range "<<(Int_t)data[bin].fTime<<ENDLOG;
158 break;
159 }
160
46fbeb8a 161 //Get the current ADC-value
162 UInt_t charge = data[bin].fCharge;
99a6d5c9 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
46fbeb8a 179 //Sum the total charge of this sequence
180 seq_charge += charge;
46fbeb8a 181 seq_average += data[bin].fTime*charge;
c3dd27a3 182 seq_error += data[bin].fTime*data[bin].fTime*charge;
183
46fbeb8a 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
99a6d5c9 198 {
adb0e845 199 LOG(AliL3Log::kFatal,"AliL3ClustFinderNew::ProcessRow","Data")
200 <<"Error in data given to the cluster finder"<<ENDLOG;
99a6d5c9 201 seq_mean = 1;
202 seq_charge = 1;
203 }
46fbeb8a 204
205 //Calculate mean in pad direction:
206 Int_t pad_mean = seq_charge*data[bin].fPad;
c3dd27a3 207 Int_t pad_error = data[bin].fPad*pad_mean;
46fbeb8a 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
c3dd27a3 216 if(difference <= fMatch) //There is a match here!!
46fbeb8a 217 {
99a6d5c9 218
219 ClusterData *local = previousPt[p];
220 if(fDeconvPad)
221 {
222 if(seq_charge > local->fLastCharge)
223 {
c3dd27a3 224 if(local->fChargeFalling) //The previous pad was falling
99a6d5c9 225 {
c3dd27a3 226 break; //create a new cluster
99a6d5c9 227 }
228 }
229 else
230 local->fChargeFalling = 1;
231 local->fLastCharge = seq_charge;
232 }
233
46fbeb8a 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:
99a6d5c9 238 //currentPt[n_current] = previousPt[p];
239
240 local->fTotalCharge += seq_charge;
241 local->fPad += pad_mean;
c3dd27a3 242 local->fPad2 += pad_error;
99a6d5c9 243 local->fTime += seq_average;
c3dd27a3 244 local->fTime2 += seq_error;
99a6d5c9 245 local->fMean = seq_mean;
246 local->fFlags++; //means we have more than one pad
247
248 currentPt[n_current] = local;
46fbeb8a 249 n_current++;
250
46fbeb8a 251 break;
c3dd27a3 252 } //Checking for match at previous pad
253 } //Loop over results on previous pad.
46fbeb8a 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.
99a6d5c9 259 //current pad will be previous pad on next pad.
46fbeb8a 260
261 //Add to the clusterlist:
262 ClusterData *tmp = &clusterlist[n_total];
263 tmp->fTotalCharge = seq_charge;
264 tmp->fPad = pad_mean;
c3dd27a3 265 tmp->fPad2 = pad_error;
46fbeb8a 266 tmp->fTime = seq_average;
c3dd27a3 267 tmp->fTime2 = seq_error;
46fbeb8a 268 tmp->fMean = seq_mean;
269 tmp->fFlags = 0; //flags for 1 pad clusters
99a6d5c9 270 if(fDeconvPad)
271 {
272 tmp->fChargeFalling = 0;
273 tmp->fLastCharge = seq_charge;
274 }
275
46fbeb8a 276 //Update list of pointers to previous pad:
277 currentPt[n_current] = &clusterlist[n_total];
278 n_total++;
279 n_current++;
280 }
99a6d5c9 281 if(fDeconvTime)
282 if(new_bin >= 0) goto redo;
283
46fbeb8a 284 }//Loop over digits on this padrow
285
286 WriteClusters(n_total,clusterlist);
46fbeb8a 287}
288
289void 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
d3b44290 298
299 Float_t xyz[3];
c3dd27a3 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;
eeddc64d 304
c3dd27a3 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
77f08cd3 312 if(fstdout==kTRUE)
c3dd27a3 313 cout<<"WriteCluster: padrow "<<fCurrentRow<<" pad "<<fpad << "+-"<<fpad2<<" time "<<ftime<<"+-"<<ftime2<<" charge "<<list[j].fTotalCharge<<endl;
d3b44290 314
494fad94 315 AliL3Transform::Slice2Sector(fCurrentSlice,fCurrentRow,thissector,thisrow);
316 AliL3Transform::Raw2Local(xyz,thissector,thisrow,fpad,ftime);
46fbeb8a 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;
d3b44290 324 }
95a00d93 325 fSpacePointData[counter].fCharge = list[j].fTotalCharge;
46fbeb8a 326 fSpacePointData[counter].fX = xyz[0];
327 fSpacePointData[counter].fY = xyz[1];
328 fSpacePointData[counter].fZ = xyz[2];
329 fSpacePointData[counter].fPadRow = fCurrentRow;
c3dd27a3 330 fSpacePointData[counter].fXYErr = fpad2;
331 fSpacePointData[counter].fZErr = ftime2;
46fbeb8a 332 fSpacePointData[counter].fID = counter
d3b44290 333 +((fCurrentSlice&0x7f)<<25)+((fCurrentPatch&0x7)<<22);//Uli
10f815d9 334#ifdef do_mc
335 Int_t trackID[3];
336 GetTrackID((Int_t)rint(fpad),(Int_t)rint(ftime),trackID);
6ab8a12a 337 //cout<<"padrow "<<fCurrentRow<<" pad "<<(Int_t)rint(fpad)<<" time "<<(Int_t)rint(ftime)<<" Trackid "<<trackID[0]<<endl;
10f815d9 338 fSpacePointData[counter].fTrackID[0] = trackID[0];
339 fSpacePointData[counter].fTrackID[1] = trackID[1];
340 fSpacePointData[counter].fTrackID[2] = trackID[2];
341#endif
46fbeb8a 342 fNClusters++;
343 counter++;
46fbeb8a 344 }
46fbeb8a 345}
10f815d9 346
347#ifdef do_mc
348void AliL3ClustFinderNew::GetTrackID(Int_t pad,Int_t time,Int_t *trackID)
349{
350 AliL3DigitRowData *rowPt = (AliL3DigitRowData*)fDigitRowData;
351
355debe1 352 trackID[0]=trackID[1]=trackID[2]=-2;
10f815d9 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;
355debe1 366 if(cpad != pad) continue;
367 if(ctime != time) continue;
368 //if(cpad != pad && ctime != ctime) continue;
10f815d9 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];
355debe1 373 break;
10f815d9 374 //cout<<"Reading trackID "<<trackID[0]<<endl;
375 }
376 break;
377 }
378
379}
380#endif