fix missing initialization of reconstructor for tracks in HLT (Kostya)
[u/mrichter/AliRoot.git] / HLT / TPCLib / AliHLTTPCClusterFinder.cxx
CommitLineData
a38a7850 1// @(#) $Id$
4aa41877 2// Original: AliHLTClustFinderNew.cxx,v 1.29 2005/06/14 10:55:21 cvetan Exp
a38a7850 3
297174de 4//**************************************************************************
5//* This file is property of and copyright by the ALICE HLT Project *
6//* ALICE Experiment at CERN, All rights reserved. *
7//* *
8//* Primary Authors: Anders Vestbo, Constantin Loizides *
9//* Developers: Kenneth Aamodt <kenneth.aamodt@student.uib.no> *
10//* Kalliopi Kanaki *
11//* for The ALICE HLT Project. *
12//* *
13//* Permission to use, copy, modify and distribute this software and its *
14//* documentation strictly for non-commercial purposes is hereby granted *
15//* without fee, provided that the above copyright notice appears in all *
16//* copies and that both the copyright notice and this permission notice *
17//* appear in the supporting documentation. The authors make no claims *
18//* about the suitability of this software for any purpose. It is *
19//* provided "as is" without express or implied warranty. *
20//**************************************************************************
a38a7850 21
600e6a1b 22/** @file AliHLTTPCClusterFinder.cxx
98034d0d 23 @author Kenneth Aamodt, Kalliopi Kanaki
600e6a1b 24 @date
25 @brief Cluster Finder for the TPC
26*/
a38a7850 27
28#include "AliHLTTPCDigitReader.h"
a38a7850 29#include "AliHLTTPCRootTypes.h"
30#include "AliHLTTPCLogging.h"
31#include "AliHLTTPCClusterFinder.h"
32#include "AliHLTTPCDigitData.h"
a38a7850 33#include "AliHLTTPCSpacePointData.h"
34#include "AliHLTTPCMemHandler.h"
46b33a24 35#include "AliHLTTPCPad.h"
8252a538 36#include <sys/time.h>
a38a7850 37
38#if __GNUC__ >= 3
39using namespace std;
40#endif
41
a38a7850 42ClassImp(AliHLTTPCClusterFinder)
43
44AliHLTTPCClusterFinder::AliHLTTPCClusterFinder()
46b33a24 45 :
2fdb1ae7 46 fClustersHWAddressVector(),
47 fRowPadVector(),
2a083ac4 48 fSpacePointData(NULL),
49 fDigitReader(NULL),
50 fPtr(NULL),
51 fSize(0),
6b15c309 52 fDeconvTime(kFALSE),
53 fDeconvPad(kFALSE),
46b33a24 54 fStdout(kFALSE),
55 fCalcerr(kTRUE),
56 fRawSP(kFALSE),
57 fFirstRow(0),
58 fLastRow(0),
2a083ac4 59 fCurrentRow(0),
60 fCurrentSlice(0),
61 fCurrentPatch(0),
62 fMatch(1),
63 fThreshold(10),
2a083ac4 64 fNClusters(0),
65 fMaxNClusters(0),
66 fXYErr(0.2),
67 fZErr(0.3),
01f43166 68 fOccupancyLimit(1.0),
64defa03 69 fUnsorted(0),
a74855c2 70 fVectorInitialized(kFALSE),
8252a538 71 fClusters(),
a74855c2 72 fNumberOfPadsInRow(NULL),
73 fNumberOfRows(0),
2fdb1ae7 74 fRowOfFirstCandidate(0),
6b15c309 75 fDoPadSelection(kFALSE),
76 fFirstTimeBin(0),
77 fLastTimeBin(AliHLTTPCTransform::GetNTimeBins()),
78 fTotalChargeOfPreviousClusterCandidate(0),
79 fChargeOfCandidatesFalling(kFALSE)
a38a7850 80{
8252a538 81 //constructor
46b33a24 82}
a38a7850 83
a38a7850 84AliHLTTPCClusterFinder::~AliHLTTPCClusterFinder()
85{
86 //destructor
8252a538 87 if(fVectorInitialized){
88 DeInitializePadArray();
89 }
90 if(fNumberOfPadsInRow){
91 delete [] fNumberOfPadsInRow;
92 fNumberOfPadsInRow=NULL;
93 }
a38a7850 94}
95
96void AliHLTTPCClusterFinder::InitSlice(Int_t slice,Int_t patch,Int_t firstrow, Int_t lastrow,Int_t nmaxpoints)
97{
98 //init slice
99 fNClusters = 0;
100 fMaxNClusters = nmaxpoints;
101 fCurrentSlice = slice;
102 fCurrentPatch = patch;
103 fFirstRow = firstrow;
104 fLastRow = lastrow;
105}
106
8252a538 107void AliHLTTPCClusterFinder::InitializePadArray(){
108 // see header file for class documentation
109
110 if(fCurrentPatch>5||fCurrentPatch<0){
111 HLTFatal("Patch is not set");
112 return;
113 }
114
115 HLTDebug("Patch number=%d",fCurrentPatch);
116
117 fFirstRow = AliHLTTPCTransform::GetFirstRow(fCurrentPatch);
118 fLastRow = AliHLTTPCTransform::GetLastRow(fCurrentPatch);
119
120 fNumberOfRows=fLastRow-fFirstRow+1;
121 fNumberOfPadsInRow= new UInt_t[fNumberOfRows];
122
123 memset( fNumberOfPadsInRow, 0, sizeof(Int_t)*(fNumberOfRows));
124
a74855c2 125 for(UInt_t i=0;i<fNumberOfRows;i++){
8252a538 126 fNumberOfPadsInRow[i]=AliHLTTPCTransform::GetNPads(i+fFirstRow);
127 AliHLTTPCPadVector tmpRow;
288f1d48 128 for(UInt_t j=0;j<=fNumberOfPadsInRow[i];j++){
8252a538 129 AliHLTTPCPad *tmpPad = new AliHLTTPCPad(2);
130 tmpPad->SetID(i,j);
131 tmpRow.push_back(tmpPad);
132 }
133 fRowPadVector.push_back(tmpRow);
134 }
135 fVectorInitialized=kTRUE;
136}
137
138Int_t AliHLTTPCClusterFinder::DeInitializePadArray()
139{
140 // see header file for class documentation
a74855c2 141 for(UInt_t i=0;i<fNumberOfRows;i++){
288f1d48 142 for(UInt_t j=0;j<=fNumberOfPadsInRow[i];j++){
8252a538 143 delete fRowPadVector[i][j];
144 fRowPadVector[i][j]=NULL;
145 }
146 fRowPadVector[i].clear();
147 }
148 fRowPadVector.clear();
149 return 1;
150}
151
152
a38a7850 153void AliHLTTPCClusterFinder::InitSlice(Int_t slice,Int_t patch,Int_t nmaxpoints)
154{
155 //init slice
156 fNClusters = 0;
157 fMaxNClusters = nmaxpoints;
158 fCurrentSlice = slice;
159 fCurrentPatch = patch;
160 fFirstRow=AliHLTTPCTransform::GetFirstRow(patch);
161 fLastRow=AliHLTTPCTransform::GetLastRow(patch);
162}
163
164void AliHLTTPCClusterFinder::SetOutputArray(AliHLTTPCSpacePointData *pt)
165{
166 //set pointer to output
167 fSpacePointData = pt;
168}
169
170void AliHLTTPCClusterFinder::Read(void* ptr,unsigned long size){
171 //set input pointer
172 fPtr = (UChar_t*)ptr;
173 fSize = size;
174}
175
176void AliHLTTPCClusterFinder::ProcessDigits()
177{
6b0fdc0e 178 int iResult=0;
a38a7850 179 bool readValue = true;
180 Int_t newRow = 0;
181 Int_t rowOffset = 0;
46b33a24 182 UShort_t time=0,newTime=0;
183 UInt_t pad=0,newPad=0;
184 AliHLTTPCSignal_t charge=0;
a38a7850 185
186 fNClusters = 0;
187
188 // initialize block for reading packed data
67fada6b 189 iResult=fDigitReader->InitBlock(fPtr,fSize,fFirstRow,fLastRow,fCurrentPatch,fCurrentSlice);
190 if (iResult<0) return;
191
a38a7850 192 readValue = fDigitReader->Next();
193
5235c3e9 194 // Matthias 08.11.2006 the following return would cause termination without writing the
6c1a9d9e 195 // ClusterData and thus would block the component. I just want to have the commented line
5235c3e9 196 // here for information
197 //if (!readValue)return;
db16520a 198
a38a7850 199 pad = fDigitReader->GetPad();
200 time = fDigitReader->GetTime();
201 fCurrentRow = fDigitReader->GetRow();
202
203 if ( fCurrentPatch >= 2 ) // Outer sector, patches 2, 3, 4, 5
204 rowOffset = AliHLTTPCTransform::GetFirstRow( 2 );
205
206 fCurrentRow += rowOffset;
207
208 UInt_t lastpad = 123456789;
2a083ac4 209 const UInt_t kPadArraySize=5000;
210 const UInt_t kClusterListSize=10000;
84645eb0 211 AliClusterData *pad1[kPadArraySize]; //2 lists for internal memory=2pads
212 AliClusterData *pad2[kPadArraySize]; //2 lists for internal memory=2pads
213 AliClusterData clusterlist[kClusterListSize]; //Clusterlist
a38a7850 214
215 AliClusterData **currentPt; //List of pointers to the current pad
216 AliClusterData **previousPt; //List of pointers to the previous pad
217 currentPt = pad2;
218 previousPt = pad1;
219 UInt_t nprevious=0,ncurrent=0,ntotal=0;
220
46b33a24 221 /* quick implementation of baseline calculation and zero suppression
222 open a pad object for each pad and delete it after processing.
223 later a list of pad objects with base line history can be used
224 The whole thing only works if we really get unprocessed raw data, if
225 the data is already zero suppressed, there might be gaps in the time
226 bins.
227 */
228 Int_t gatingGridOffset=50;
6b15c309 229 if(fFirstTimeBin>0){
230 gatingGridOffset=fFirstTimeBin;
231 }
46b33a24 232 AliHLTTPCPad baseline(gatingGridOffset, AliHLTTPCTransform::GetNTimeBins());
233 // just to make later conversion to a list of objects easier
84645eb0 234 AliHLTTPCPad* pCurrentPad=NULL;
6b15c309 235 /*
236 if (fSignalThreshold>=0) {
84645eb0 237 pCurrentPad=&baseline;
238 baseline.SetThreshold(fSignalThreshold);
239 }
6b15c309 240 */
6b0fdc0e 241 while ( readValue!=0 && iResult>=0){ // Reads through all digits in block
242 iResult=0;
a38a7850 243
244 if(pad != lastpad){
245 //This is a new pad
246
247 //Switch the lists:
248 if(currentPt == pad2){
249 currentPt = pad1;
250 previousPt = pad2;
251 }
252 else {
253 currentPt = pad2;
254 previousPt = pad1;
255 }
256 nprevious = ncurrent;
257 ncurrent = 0;
258 if(pad != lastpad+1){
259 //this happens if there is a pad with no signal.
260 nprevious = ncurrent = 0;
261 }
262 lastpad = pad;
263 }
264
265 Bool_t newcluster = kTRUE;
266 UInt_t seqcharge=0,seqaverage=0,seqerror=0;
2a083ac4 267 AliHLTTPCSignal_t lastcharge=0;
268 UInt_t bLastWasFalling=0;
a38a7850 269 Int_t newbin=-1;
270
271
272 if(fDeconvTime){
273 redo: //This is a goto.
274
275 if(newbin > -1){
276 //bin = newbin;
277 newbin = -1;
278 }
279
280 lastcharge=0;
2a083ac4 281 bLastWasFalling = 0;
a38a7850 282 }
283
6b0fdc0e 284 while(iResult>=0){ //Loop over time bins of current pad
46b33a24 285 // read all the values for one pad at once to calculate the base line
286 if (pCurrentPad) {
287 if (!pCurrentPad->IsStarted()) {
288 //HLTDebug("reading data for pad %d, padrow %d", fDigitReader->GetPad(), fDigitReader->GetRow()+rowOffset);
289 pCurrentPad->SetID(fDigitReader->GetRow()+rowOffset,fDigitReader->GetPad());
290 if ((pCurrentPad->StartEvent())>=0) {
291 do {
292 if ((fDigitReader->GetRow()+rowOffset)!=pCurrentPad->GetRowNumber()) break;
293 if (fDigitReader->GetPad()!=pCurrentPad->GetPadNumber()) break;
294 pCurrentPad->SetRawData(fDigitReader->GetTime(), fDigitReader->GetSignal());
295 //HLTDebug("set raw data to pad: bin %d charge %d", fDigitReader->GetTime(), fDigitReader->GetSignal());
296 } while ((readValue = fDigitReader->Next())!=0);
297 }
298 pCurrentPad->CalculateBaseLine(AliHLTTPCTransform::GetNTimeBins()/2);
299 if (pCurrentPad->Next(kTRUE/*do zero suppression*/)==0) {
84645eb0 300 //HLTDebug("no data available after zero suppression");
46b33a24 301 pCurrentPad->StopEvent();
302 pCurrentPad->ResetHistory();
303 break;
304 }
305 time=pCurrentPad->GetCurrentPosition();
306 if (time>pCurrentPad->GetSize()) {
307 HLTError("invalid time bin for pad");
308 break;
309 }
310 }
311 }
46b33a24 312 if (pCurrentPad) {
5235c3e9 313 Float_t occupancy=pCurrentPad->GetOccupancy();
314 //HLTDebug("pad %d occupancy level: %f", pCurrentPad->GetPadNumber(), occupancy);
315 if ( occupancy < fOccupancyLimit ) {
316 charge = pCurrentPad->GetCorrectedData();
317 } else {
318 charge = 0;
319 //HLTDebug("ignoring pad %d with occupancy level %f", pCurrentPad->GetPadNumber(), occupancy);
320 }
46b33a24 321 } else {
5235c3e9 322 charge = fDigitReader->GetSignal();
46b33a24 323 }
324 //HLTDebug("get next charge value: position %d charge %d", time, charge);
a38a7850 325
738c049f 326
327 // CHARGE DEBUG
328 if (fDigitReader->GetRow() == 90){
329///// LOG(AliHLTTPCLog::kFatal,"AliHLTTPCClusterFinder::Row","row90") << "PAD=" << fDigitReader->GetPad() << " TIME=" << fDigitReader->GetTime()
330 // << " SIGNAL=" << fDigitReader->GetSignal() << ENDLOG;
331
332 }
333
a38a7850 334 if(time >= AliHLTTPCTransform::GetNTimeBins()){
6b0fdc0e 335 HLTWarning("Pad %d: Timebin (%d) out of range (%d)", pad, time, AliHLTTPCTransform::GetNTimeBins());
336 iResult=-ERANGE;
a38a7850 337 }
46b33a24 338
339
a38a7850 340 //Get the current ADC-value
341 if(fDeconvTime){
342
343 //Check if the last pixel in the sequence is smaller than this
344 if(charge > lastcharge){
2a083ac4 345 if(bLastWasFalling){
a38a7850 346 newbin = 1;
347 break;
348 }
349 }
2a083ac4 350 else bLastWasFalling = 1; //last pixel was larger than this
a38a7850 351 lastcharge = charge;
352 }
353
354 //Sum the total charge of this sequence
355 seqcharge += charge;
356 seqaverage += time*charge;
357 seqerror += time*time*charge;
358
46b33a24 359 if (pCurrentPad) {
360
361 if((pCurrentPad->Next(kTRUE/*do zero suppression*/))==0) {
362 pCurrentPad->StopEvent();
363 pCurrentPad->ResetHistory();
364 if(readValue) {
365 newPad = fDigitReader->GetPad();
366 newTime = fDigitReader->GetTime();
367 newRow = fDigitReader->GetRow() + rowOffset;
368 }
369 break;
370 }
371
372 newPad=pCurrentPad->GetPadNumber();
373 newTime=pCurrentPad->GetCurrentPosition();
374 newRow=pCurrentPad->GetRowNumber();
375 } else {
a38a7850 376 readValue = fDigitReader->Next();
a38a7850 377 //Check where to stop:
378 if(!readValue) break; //No more value
379
380 newPad = fDigitReader->GetPad();
381 newTime = fDigitReader->GetTime();
382 newRow = fDigitReader->GetRow() + rowOffset;
46b33a24 383 }
a38a7850 384
385 if(newPad != pad)break; //new pad
386 if(newTime != time+1) break; //end of sequence
6b0fdc0e 387 if(iResult<0) break;
a38a7850 388
389 // pad = newpad; is equal
390 time = newTime;
391
392 }//end loop over sequence
393
46b33a24 394 //HLTDebug("ended time bin sequence loop: seqcharge=%d readValue=%d", seqcharge, readValue);
395 //HLTDebug("pad=%d newpad=%d current row=%d newrow=%d", pad, newPad, fCurrentRow, newRow);
396 if (seqcharge<=0) {
397 // with active zero suppression zero values are possible
398 continue;
399 }
a38a7850 400
401 //Calculate mean of sequence:
402 Int_t seqmean=0;
403 if(seqcharge)
404 seqmean = seqaverage/seqcharge;
405 else{
406 LOG(AliHLTTPCLog::kFatal,"AliHLTTPCClusterFinder::ProcessRow","Data")
407 <<"Error in data given to the cluster finder"<<ENDLOG;
408 seqmean = 1;
409 seqcharge = 1;
410 }
411
412 //Calculate mean in pad direction:
413 Int_t padmean = seqcharge*pad;
414 Int_t paderror = pad*padmean;
415
a38a7850 416 //Compare with results on previous pad:
84645eb0 417 for(UInt_t p=0; p<nprevious && p<kPadArraySize && ncurrent<kPadArraySize; p++){
a38a7850 418
419 //dont merge sequences on the same pad twice
420 if(previousPt[p]->fLastMergedPad==pad) continue;
421
422 Int_t difference = seqmean - previousPt[p]->fMean;
423 if(difference < -fMatch) break;
424
425 if(difference <= fMatch){ //There is a match here!!
426 AliClusterData *local = previousPt[p];
427
428 if(fDeconvPad){
429 if(seqcharge > local->fLastCharge){
430 if(local->fChargeFalling){ //The previous pad was falling
431 break; //create a new cluster
432 }
433 }
434 else local->fChargeFalling = 1;
435 local->fLastCharge = seqcharge;
436 }
437
438 //Don't create a new cluster, because we found a match
439 newcluster = kFALSE;
440
441 //Update cluster on current pad with the matching one:
442 local->fTotalCharge += seqcharge;
443 local->fPad += padmean;
444 local->fPad2 += paderror;
445 local->fTime += seqaverage;
446 local->fTime2 += seqerror;
447 local->fMean = seqmean;
448 local->fFlags++; //means we have more than one pad
449 local->fLastMergedPad = pad;
450
451 currentPt[ncurrent] = local;
452 ncurrent++;
453
454 break;
455 } //Checking for match at previous pad
456 } //Loop over results on previous pad.
457
84645eb0 458 if(newcluster && ncurrent<kPadArraySize){
a38a7850 459 //Start a new cluster. Add it to the clusterlist, and update
460 //the list of pointers to clusters in current pad.
461 //current pad will be previous pad on next pad.
462
463 //Add to the clusterlist:
464 AliClusterData *tmp = &clusterlist[ntotal];
465 tmp->fTotalCharge = seqcharge;
466 tmp->fPad = padmean;
467 tmp->fPad2 = paderror;
468 tmp->fTime = seqaverage;
469 tmp->fTime2 = seqerror;
470 tmp->fMean = seqmean;
471 tmp->fFlags = 0; //flags for single pad clusters
472 tmp->fLastMergedPad = pad;
473
474 if(fDeconvPad){
475 tmp->fChargeFalling = 0;
476 tmp->fLastCharge = seqcharge;
477 }
478
479 //Update list of pointers to previous pad:
480 currentPt[ncurrent] = &clusterlist[ntotal];
481 ntotal++;
482 ncurrent++;
483 }
484
485 if(fDeconvTime)
486 if(newbin >= 0) goto redo;
487
488 // to prevent endless loop
489 if(time >= AliHLTTPCTransform::GetNTimeBins()){
46b33a24 490 HLTWarning("Timebin (%d) out of range (%d)", time, AliHLTTPCTransform::GetNTimeBins());
a38a7850 491 break;
492 }
493
494
495 if(!readValue) break; //No more value
84645eb0 496
497 if (ntotal>=kClusterListSize || ncurrent>=kPadArraySize) {
498 HLTWarning("pad array size exceeded ntotal=%d ncurrent=%d, skip rest of the data", ntotal, ncurrent);
499 break;
500 }
a38a7850 501
502 if(fCurrentRow != newRow){
503 WriteClusters(ntotal,clusterlist);
504
505 lastpad = 123456789;
506
507 currentPt = pad2;
508 previousPt = pad1;
509 nprevious=0;
510 ncurrent=0;
511 ntotal=0;
512
513 fCurrentRow = newRow;
514 }
515
516 pad = newPad;
517 time = newTime;
db16520a 518
a38a7850 519 } // END while(readValue)
520
521 WriteClusters(ntotal,clusterlist);
522
46b33a24 523 HLTInfo("ClusterFinder found %d clusters in slice %d patch %d", fNClusters, fCurrentSlice, fCurrentPatch);
a38a7850 524
525} // ENDEND
526
527void AliHLTTPCClusterFinder::WriteClusters(Int_t nclusters,AliClusterData *list)
528{
529 //write cluster to output pointer
8f8bf0af 530 Int_t thisrow=-1,thissector=-1;
a38a7850 531 UInt_t counter = fNClusters;
532
533 for(int j=0; j<nclusters; j++)
534 {
d8b2f74a 535
536
537
a38a7850 538 if(!list[j].fFlags) continue; //discard single pad clusters
539 if(list[j].fTotalCharge < fThreshold) continue; //noise cluster
540
541 Float_t xyz[3];
542 Float_t fpad =(Float_t)list[j].fPad / list[j].fTotalCharge;
543 Float_t fpad2=fXYErr*fXYErr; //fixed given error
544 Float_t ftime =(Float_t)list[j].fTime / list[j].fTotalCharge;
545 Float_t ftime2=fZErr*fZErr; //fixed given error
546
738c049f 547
aff6e981 548
549 if(fUnsorted){
550 fCurrentRow=list[j].fRow;
551 }
552
738c049f 553
a38a7850 554 if(fCalcerr) { //calc the errors, otherwice take the fixed error
555 Int_t patch = AliHLTTPCTransform::GetPatch(fCurrentRow);
556 UInt_t q2=list[j].fTotalCharge*list[j].fTotalCharge;
5611e981 557 // Float_t sy2=list[j].fPad2 * list[j].fTotalCharge - list[j].fPad * list[j].fPad;
558 Float_t sy2=(Float_t)list[j].fPad2 * list[j].fTotalCharge - (Float_t)list[j].fPad * list[j].fPad;
a38a7850 559 sy2/=q2;
560 if(sy2 < 0) {
561 LOG(AliHLTTPCLog::kError,"AliHLTTPCClusterFinder::WriteClusters","Cluster width")
562 <<"SigmaY2 negative "<<sy2<<" on row "<<fCurrentRow<<" "<<fpad<<" "<<ftime<<ENDLOG;
563 continue;
564 } else {
565 if(!fRawSP){
566 fpad2 = (sy2 + 1./12)*AliHLTTPCTransform::GetPadPitchWidth(patch)*AliHLTTPCTransform::GetPadPitchWidth(patch);
567 if(sy2 != 0){
568 fpad2*=0.108; //constants are from offline studies
569 if(patch<2)
570 fpad2*=2.07;
571 }
572 } else fpad2=sy2; //take the width not the error
573 }
5611e981 574 // Float_t sz2=list[j].fTime2*list[j].fTotalCharge - list[j].fTime*list[j].fTime;
575 Float_t sz2=(Float_t)list[j].fTime2*list[j].fTotalCharge - (Float_t)list[j].fTime*list[j].fTime;
a38a7850 576 sz2/=q2;
577 if(sz2 < 0){
578 LOG(AliHLTTPCLog::kError,"AliHLTTPCClusterFinder::WriteClusters","Cluster width")
579 <<"SigmaZ2 negative "<<sz2<<" on row "<<fCurrentRow<<" "<<fpad<<" "<<ftime<<ENDLOG;
580 continue;
581 } else {
582 if(!fRawSP){
583 ftime2 = (sz2 + 1./12)*AliHLTTPCTransform::GetZWidth()*AliHLTTPCTransform::GetZWidth();
584 if(sz2 != 0) {
585 ftime2 *= 0.169; //constants are from offline studies
586 if(patch<2)
587 ftime2 *= 1.77;
588 }
589 } else ftime2=sz2; //take the width, not the error
590 }
591 }
592 if(fStdout==kTRUE)
8252a538 593 HLTInfo("WriteCluster: padrow %d pad %d +- %d time +- %d charge %d",fCurrentRow, fpad, fpad2, ftime, ftime2, list[j].fTotalCharge);
a38a7850 594
595 if(!fRawSP){
596 AliHLTTPCTransform::Slice2Sector(fCurrentSlice,fCurrentRow,thissector,thisrow);
597 AliHLTTPCTransform::Raw2Local(xyz,thissector,thisrow,fpad,ftime);
598
599 if(xyz[0]==0) LOG(AliHLTTPCLog::kError,"AliHLTTPCClustFinder","Cluster Finder")
600 <<AliHLTTPCLog::kDec<<"Zero cluster"<<ENDLOG;
601 if(fNClusters >= fMaxNClusters)
602 {
603 LOG(AliHLTTPCLog::kError,"AliHLTTPCClustFinder::WriteClusters","Cluster Finder")
604 <<AliHLTTPCLog::kDec<<"Too many clusters "<<fNClusters<<ENDLOG;
605 return;
606 }
607
608 fSpacePointData[counter].fX = xyz[0];
6b15c309 609 // fSpacePointData[counter].fY = xyz[1];
610 if(fCurrentSlice<18){
611 fSpacePointData[counter].fY = xyz[1];
612 }
613 else{
614 fSpacePointData[counter].fY = -1*xyz[1];
615 }
a38a7850 616 fSpacePointData[counter].fZ = xyz[2];
6b15c309 617
a38a7850 618 } else {
619 fSpacePointData[counter].fX = fCurrentRow;
620 fSpacePointData[counter].fY = fpad;
621 fSpacePointData[counter].fZ = ftime;
622 }
623
624 fSpacePointData[counter].fCharge = list[j].fTotalCharge;
625 fSpacePointData[counter].fPadRow = fCurrentRow;
626 fSpacePointData[counter].fSigmaY2 = fpad2;
627 fSpacePointData[counter].fSigmaZ2 = ftime2;
628
6b15c309 629 fSpacePointData[counter].fMaxQ = list[j].fQMax;
630
44be0fde 631 fSpacePointData[counter].fUsed = kFALSE; // only used / set in AliHLTTPCDisplay
738c049f 632 fSpacePointData[counter].fTrackN = -1; // only used / set in AliHLTTPCDisplay
44be0fde 633
a38a7850 634 Int_t patch=fCurrentPatch;
635 if(patch==-1) patch=0; //never store negative patch number
636 fSpacePointData[counter].fID = counter
637 +((fCurrentSlice&0x7f)<<25)+((patch&0x7)<<22);//Uli
638
639#ifdef do_mc
640 Int_t trackID[3];
641 GetTrackID((Int_t)rint(fpad),(Int_t)rint(ftime),trackID);
642
643 fSpacePointData[counter].fTrackID[0] = trackID[0];
644 fSpacePointData[counter].fTrackID[1] = trackID[1];
645 fSpacePointData[counter].fTrackID[2] = trackID[2];
646
a38a7850 647#endif
648
649 fNClusters++;
650 counter++;
651 }
652}
653
654// STILL TO FIX ----------------------------------------------------------------------------
655
656#ifdef do_mc
657void AliHLTTPCClusterFinder::GetTrackID(Int_t pad,Int_t time,Int_t *trackID)
658{
659 //get mc id
660 AliHLTTPCDigitRowData *rowPt = (AliHLTTPCDigitRowData*)fDigitRowData;
661
662 trackID[0]=trackID[1]=trackID[2]=-2;
a38a7850 663 for(Int_t i=fFirstRow; i<=fLastRow; i++){
664 if(rowPt->fRow < (UInt_t)fCurrentRow){
665 AliHLTTPCMemHandler::UpdateRowPointer(rowPt);
666 continue;
667 }
668 AliHLTTPCDigitData *digPt = (AliHLTTPCDigitData*)rowPt->fDigitData;
669 for(UInt_t j=0; j<rowPt->fNDigit; j++){
670 Int_t cpad = digPt[j].fPad;
671 Int_t ctime = digPt[j].fTime;
672 if(cpad != pad) continue;
673 if(ctime != time) continue;
674
675 trackID[0] = digPt[j].fTrackID[0];
676 trackID[1] = digPt[j].fTrackID[1];
677 trackID[2] = digPt[j].fTrackID[2];
678
a38a7850 679 break;
680 }
681 break;
682 }
683}
684#endif
01f43166 685
686//----------------------------------Methods for the new unsorted way of reading the data --------------------------------
687
a74855c2 688void AliHLTTPCClusterFinder::ReadDataUnsorted(void* ptr,unsigned long size)
64defa03 689{
01f43166 690 //set input pointer
691 fPtr = (UChar_t*)ptr;
692 fSize = size;
693
8252a538 694 if(!fVectorInitialized){
695 InitializePadArray();
696 }
b1c46961 697
5863c71a 698 if (fDigitReader->InitBlock(fPtr,fSize,fFirstRow,fLastRow,fCurrentPatch,fCurrentSlice)<0) {
699 HLTError("failed setting up digit reader (InitBlock)");
700 return;
701 }
f5b70376 702
8252a538 703 while(fDigitReader->NextChannel()){
704 UInt_t row=fDigitReader->GetRow();
705 UInt_t pad=fDigitReader->GetPad();
70a4e20c 706
3789100d 707 if(row>=fRowPadVector.size()){
708 HLTError("Row number is to large: %d, max is %d",row,fRowPadVector.size()-1);
709 continue;
710 }
711 if(pad>=fRowPadVector[row].size()){
712 HLTError("Pad number is to large: %d, max is %d",pad,fRowPadVector[row].size());
713 continue;
714 }
715
8252a538 716 while(fDigitReader->NextBunch()){
717 if(fDigitReader->GetBunchSize()>1){//to remove single timebin values, this will have to change at some point
8252a538 718 UInt_t time = fDigitReader->GetTime();
6b15c309 719 if((Int_t)time>=fFirstTimeBin && (Int_t)time+fDigitReader->GetBunchSize()<=fLastTimeBin){
720 const UInt_t *bunchData= fDigitReader->GetSignals();
721 AliHLTTPCClusters candidate;
722 for(Int_t i=0;i<fDigitReader->GetBunchSize();i++){
723 candidate.fTotalCharge+=bunchData[i];
724 candidate.fTime += time*bunchData[i];
725 candidate.fTime2 += time*time*bunchData[i];
726 if(bunchData[i]>candidate.fQMax){
727 candidate.fQMax=bunchData[i];
728 }
729 time++;
730 }
731 if(candidate.fTotalCharge>0){
732 candidate.fMean=candidate.fTime/candidate.fTotalCharge;
733 candidate.fPad=candidate.fTotalCharge*pad;
734 candidate.fPad2=candidate.fPad*pad;
735 candidate.fLastMergedPad=pad;
736 candidate.fRowNumber=row+fDigitReader->GetRowOffset();
737 }
3789100d 738 if(fRowPadVector[row][pad] != NULL){
739 fRowPadVector[row][pad]->AddClusterCandidate(candidate);
740 }
8252a538 741 }
6b15c309 742 }
743 }
744 }
745}
746
747void AliHLTTPCClusterFinder::ReadDataUnsortedDeconvoluteTime(void* ptr,unsigned long size)
748{
749 //set input pointer
750 fPtr = (UChar_t*)ptr;
751 fSize = size;
752
753 if(!fVectorInitialized){
754 InitializePadArray();
755 }
756
5863c71a 757 if (fDigitReader->InitBlock(fPtr,fSize,fFirstRow,fLastRow,fCurrentPatch,fCurrentSlice)<0) {
758 HLTError("failed setting up digit reader (InitBlock)");
759 return;
760 }
6b15c309 761
762 while(fDigitReader->NextChannel()){
763 UInt_t row=fDigitReader->GetRow();
764 UInt_t pad=fDigitReader->GetPad();
765
766 while(fDigitReader->NextBunch()){
767 if(fDigitReader->GetBunchSize()>1){//to remove single timebin values, this will have to change at some point
768 UInt_t time = fDigitReader->GetTime();
769 if((Int_t)time>=fFirstTimeBin && (Int_t)time+fDigitReader->GetBunchSize()<=fLastTimeBin){
770 Int_t indexInBunchData=0;
771 Bool_t moreDataInBunch=kFALSE;
772 UInt_t prevSignal=0;
773 Bool_t signalFalling=kFALSE;
774 const UInt_t *bunchData= fDigitReader->GetSignals();
775 do{
776 AliHLTTPCClusters candidate;
777 for(Int_t i=indexInBunchData;i<fDigitReader->GetBunchSize();i++){
778 // Checks if one need to deconvolute the signals
779 if(bunchData[i]>prevSignal && signalFalling==kTRUE){
780 if(i<fDigitReader->GetBunchSize()-1){ // means there are more than one signal left in the bunch
781 moreDataInBunch=kTRUE;
782 prevSignal=0;
783 }
784 break;
785 }
786
787 // Checks if the signal is 0, then quit processing the data.
788 if(bunchData[i]==0 && i<fDigitReader->GetBunchSize()-1){//means we have 0 data fom the rcu, might happen depending on the configuration settings
789 moreDataInBunch=kTRUE;
790 prevSignal=0;
791 break;
792 }
793
794 if(prevSignal>bunchData[i]){//means the peak of the signal has been reached and deconvolution will happen if the signal rise again.
795 signalFalling=kTRUE;
796 }
797 candidate.fTotalCharge+=bunchData[i];
798 candidate.fTime += time*bunchData[i];
799 candidate.fTime2 += time*time*bunchData[i];
800 if(bunchData[i]>candidate.fQMax){
801 candidate.fQMax=bunchData[i];
802 }
803
804 prevSignal=bunchData[i];
805 time++;
806 indexInBunchData++;
807 }
808 if(candidate.fTotalCharge>0){
809 candidate.fMean=candidate.fTime/candidate.fTotalCharge;
810 candidate.fPad=candidate.fTotalCharge*pad;
811 candidate.fPad2=candidate.fPad*pad;
812 candidate.fLastMergedPad=pad;
813 candidate.fRowNumber=row+fDigitReader->GetRowOffset();
814 }
815 fRowPadVector[row][pad]->AddClusterCandidate(candidate);
816 if(indexInBunchData<fDigitReader->GetBunchSize()-1){
817 moreDataInBunch=kFALSE;
818 }
819 }while(moreDataInBunch);
8252a538 820 }
8252a538 821 }
822 }
b1c46961 823 }
8252a538 824}
825
826Bool_t AliHLTTPCClusterFinder::ComparePads(AliHLTTPCPad *nextPad,AliHLTTPCClusters* cluster,Int_t nextPadToRead){
827 //Checking if we have a match on the next pad
a74855c2 828 for(UInt_t candidateNumber=0;candidateNumber<nextPad->fClusterCandidates.size();candidateNumber++){
8252a538 829 AliHLTTPCClusters *candidate =&nextPad->fClusterCandidates[candidateNumber];
830 if(cluster->fMean-candidate->fMean==1 || candidate->fMean-cluster->fMean==1 || cluster->fMean-candidate->fMean==0){
6b15c309 831 if(fDeconvPad){
832 if(candidate->fTotalCharge<fTotalChargeOfPreviousClusterCandidate){//peak is reached
833 fChargeOfCandidatesFalling=kTRUE;
834 }
835 if(candidate->fTotalCharge>fTotalChargeOfPreviousClusterCandidate && fChargeOfCandidatesFalling==kTRUE){//we have deconvolution
836 return kFALSE;
837 }
838 }
8252a538 839 cluster->fMean=candidate->fMean;
840 cluster->fTotalCharge+=candidate->fTotalCharge;
841 cluster->fTime += candidate->fTime;
842 cluster->fTime2 += candidate->fTime2;
843 cluster->fPad+=candidate->fPad;
5611e981 844 cluster->fPad2+=candidate->fPad2;
8252a538 845 cluster->fLastMergedPad=candidate->fPad;
6b15c309 846 if(candidate->fQMax>cluster->fQMax){
847 cluster->fQMax=candidate->fQMax;
848 }
849
2fdb1ae7 850 if(fDoPadSelection){
6b15c309 851 UInt_t rowNo = nextPad->GetRowNumber();
852 UInt_t padNo = nextPad->GetPadNumber();
853 if(padNo-1>0){
854 fRowPadVector[rowNo][padNo-2]->fSelectedPad=kTRUE;
855 fRowPadVector[rowNo][padNo-2]->fHWAddress=(AliHLTUInt16_t)fDigitReader->GetAltroBlockHWaddr(rowNo,padNo-2);
856 }
857 fRowPadVector[rowNo][padNo-1]->fSelectedPad=kTRUE;// quick solution to set the first pad to selected
858 fRowPadVector[rowNo][padNo-1]->fHWAddress=(AliHLTUInt16_t)fDigitReader->GetAltroBlockHWaddr(rowNo,padNo-1);
859 fRowPadVector[rowNo][padNo]->fSelectedPad=kTRUE;
860 fRowPadVector[rowNo][padNo]->fHWAddress=(AliHLTUInt16_t)fDigitReader->GetAltroBlockHWaddr(rowNo,padNo);
2fdb1ae7 861 }
8252a538 862
863 //setting the matched pad to used
864 nextPad->fUsedClusterCandidates[candidateNumber]=1;
865 nextPadToRead++;
a74855c2 866 if(nextPadToRead<(Int_t)fNumberOfPadsInRow[fRowOfFirstCandidate]){
8252a538 867 nextPad=fRowPadVector[fRowOfFirstCandidate][nextPadToRead];
868 ComparePads(nextPad,cluster,nextPadToRead);
869 }
870 else{
871 return kFALSE;
872 }
873 }
874 else{
875 return kFALSE;
876 }
b1c46961 877 }
8252a538 878 return kFALSE;
01f43166 879}
64defa03 880
2fdb1ae7 881Int_t AliHLTTPCClusterFinder::FillHWAddressList(AliHLTUInt16_t *hwaddlist, Int_t maxHWadd){
882 Int_t counter=0;
883 for(UInt_t row=0;row<fNumberOfRows;row++){
884 for(UInt_t pad=0;pad<fNumberOfPadsInRow[row]-1;pad++){
885 if(fRowPadVector[row][pad]->fSelectedPad){
886 if(counter<maxHWadd){
887 hwaddlist[counter]=(AliHLTUInt16_t)fRowPadVector[row][pad]->fHWAddress;
2fdb1ae7 888 counter++;
889 }
890 else{
891 HLTWarning("To many hardwareaddresses, skip adding");
892 }
893
894 }
895 }
896 }
897 return counter;
898}
899
64defa03 900void AliHLTTPCClusterFinder::FindClusters()
901{
b1c46961 902 // see header file for function documentation
8252a538 903
904 AliHLTTPCClusters* tmpCandidate=NULL;
a74855c2 905 for(UInt_t row=0;row<fNumberOfRows;row++){
8252a538 906 fRowOfFirstCandidate=row;
288f1d48 907 for(UInt_t pad=0;pad<fNumberOfPadsInRow[row];pad++){
8252a538 908 AliHLTTPCPad *tmpPad=fRowPadVector[row][pad];
909 for(size_t candidate=0;candidate<tmpPad->fClusterCandidates.size();candidate++){
910 if(tmpPad->fUsedClusterCandidates[candidate]){
911 continue;
912 }
913 tmpCandidate=&tmpPad->fClusterCandidates[candidate];
914 UInt_t tmpTotalCharge=tmpCandidate->fTotalCharge;
6b15c309 915
8252a538 916 ComparePads(fRowPadVector[row][pad+1],tmpCandidate,pad+1);
917 if(tmpCandidate->fTotalCharge>tmpTotalCharge){
918 //we have a cluster
919 fClusters.push_back(*tmpCandidate);
920 }
921 }
368ff9d3 922 tmpPad->ClearCandidates();
8252a538 923 }
288f1d48 924 fRowPadVector[row][fNumberOfPadsInRow[row]]->ClearCandidates();
01f43166 925 }
8252a538 926
927 HLTInfo("Found %d clusters.",fClusters.size());
928
929 //TODO: Change so it stores AliHLTTPCSpacePointData directly, instead of this copying
930
931 AliClusterData * clusterlist = new AliClusterData[fClusters.size()]; //Clusterlist
932 for(unsigned int i=0;i<fClusters.size();i++){
933 clusterlist[i].fTotalCharge = fClusters[i].fTotalCharge;
934 clusterlist[i].fPad = fClusters[i].fPad;
935 clusterlist[i].fPad2 = fClusters[i].fPad2;
936 clusterlist[i].fTime = fClusters[i].fTime;
937 clusterlist[i].fTime2 = fClusters[i].fTime2;
938 clusterlist[i].fMean = fClusters[i].fMean;
939 clusterlist[i].fFlags = fClusters[i].fFlags;
940 clusterlist[i].fChargeFalling = fClusters[i].fChargeFalling;
941 clusterlist[i].fLastCharge = fClusters[i].fLastCharge;
942 clusterlist[i].fLastMergedPad = fClusters[i].fLastMergedPad;
943 clusterlist[i].fRow = fClusters[i].fRowNumber;
fec96a66 944 clusterlist[i].fQMax = fClusters[i].fQMax;
945
8252a538 946 }
947 // PrintClusters();
948 WriteClusters(fClusters.size(),clusterlist);
e83e889b 949 delete [] clusterlist;
8252a538 950 fClusters.clear();
01f43166 951}
64defa03 952
8252a538 953void AliHLTTPCClusterFinder::PrintClusters()
954{
955 // see header file for class documentation
956 for(size_t i=0;i<fClusters.size();i++){
957 HLTInfo("Cluster number: %d",i);
6b15c309 958 HLTInfo("Row: %d \t Pad: %d",fClusters[i].fRowNumber,fClusters[i].fPad/fClusters[i].fTotalCharge);
8252a538 959 HLTInfo("Total Charge: %d",fClusters[i].fTotalCharge);
960 HLTInfo("fPad: %d",fClusters[i].fPad);
961 HLTInfo("PadError: %d",fClusters[i].fPad2);
6b15c309 962 HLTInfo("TimeMean: %d",fClusters[i].fTime/fClusters[i].fTotalCharge);
8252a538 963 HLTInfo("TimeError: %d",fClusters[i].fTime2);
964 HLTInfo("EndOfCluster:");
965 }
966}
967
01f43166 968void AliHLTTPCClusterFinder::WriteClusters(Int_t nclusters,AliHLTTPCClusters *list)//This is used when using the AliHLTTPCClusters class for cluster data
969{
970 //write cluster to output pointer
971 Int_t thisrow,thissector;
972 UInt_t counter = fNClusters;
973
974 for(int j=0; j<nclusters; j++)
975 {
976 if(!list[j].fFlags) continue; //discard single pad clusters
977 if(list[j].fTotalCharge < fThreshold) continue; //noise cluster
978
979 Float_t xyz[3];
980 Float_t fpad =(Float_t)list[j].fPad / list[j].fTotalCharge;
981 Float_t fpad2=fXYErr*fXYErr; //fixed given error
982 Float_t ftime =(Float_t)list[j].fTime / list[j].fTotalCharge;
983 Float_t ftime2=fZErr*fZErr; //fixed given error
984
d8b2f74a 985
01f43166 986 if(fCalcerr) { //calc the errors, otherwice take the fixed error
987 Int_t patch = AliHLTTPCTransform::GetPatch(fCurrentRow);
988 UInt_t q2=list[j].fTotalCharge*list[j].fTotalCharge;
989 Float_t sy2=list[j].fPad2 * list[j].fTotalCharge - list[j].fPad * list[j].fPad;
990 sy2/=q2;
991 if(sy2 < 0) {
992 LOG(AliHLTTPCLog::kError,"AliHLTTPCClusterFinder::WriteClusters","Cluster width")
993 <<"SigmaY2 negative "<<sy2<<" on row "<<fCurrentRow<<" "<<fpad<<" "<<ftime<<ENDLOG;
994 continue;
995 } else {
996 if(!fRawSP){
997 fpad2 = (sy2 + 1./12)*AliHLTTPCTransform::GetPadPitchWidth(patch)*AliHLTTPCTransform::GetPadPitchWidth(patch);
998 if(sy2 != 0){
999 fpad2*=0.108; //constants are from offline studies
1000 if(patch<2)
1001 fpad2*=2.07;
1002 }
1003 } else fpad2=sy2; //take the width not the error
1004 }
1005 Float_t sz2=list[j].fTime2*list[j].fTotalCharge - list[j].fTime*list[j].fTime;
1006 sz2/=q2;
1007 if(sz2 < 0){
1008 LOG(AliHLTTPCLog::kError,"AliHLTTPCClusterFinder::WriteClusters","Cluster width")
1009 <<"SigmaZ2 negative "<<sz2<<" on row "<<fCurrentRow<<" "<<fpad<<" "<<ftime<<ENDLOG;
1010 continue;
1011 } else {
1012 if(!fRawSP){
1013 ftime2 = (sz2 + 1./12)*AliHLTTPCTransform::GetZWidth()*AliHLTTPCTransform::GetZWidth();
1014 if(sz2 != 0) {
1015 ftime2 *= 0.169; //constants are from offline studies
1016 if(patch<2)
1017 ftime2 *= 1.77;
1018 }
1019 } else ftime2=sz2; //take the width, not the error
1020 }
1021 }
1022 if(fStdout==kTRUE)
8252a538 1023 HLTInfo("WriteCluster: padrow %d pad %d +- %d time +- %d charge %d",fCurrentRow, fpad, fpad2, ftime, ftime2, list[j].fTotalCharge);
1024
01f43166 1025 if(!fRawSP){
1026 AliHLTTPCTransform::Slice2Sector(fCurrentSlice,fCurrentRow,thissector,thisrow);
1027 AliHLTTPCTransform::Raw2Local(xyz,thissector,thisrow,fpad,ftime);
1028
1029 if(xyz[0]==0) LOG(AliHLTTPCLog::kError,"AliHLTTPCClustFinder","Cluster Finder")
1030 <<AliHLTTPCLog::kDec<<"Zero cluster"<<ENDLOG;
1031 if(fNClusters >= fMaxNClusters)
1032 {
1033 LOG(AliHLTTPCLog::kError,"AliHLTTPCClustFinder::WriteClusters","Cluster Finder")
1034 <<AliHLTTPCLog::kDec<<"Too many clusters "<<fNClusters<<ENDLOG;
1035 return;
1036 }
1037
1038 fSpacePointData[counter].fX = xyz[0];
6b15c309 1039 // fSpacePointData[counter].fY = xyz[1];
1040 if(fCurrentSlice<18){
1041 fSpacePointData[counter].fY = xyz[1];
1042 }
1043 else{
1044 fSpacePointData[counter].fY = -1*xyz[1];
1045 }
01f43166 1046 fSpacePointData[counter].fZ = xyz[2];
1047
1048 } else {
1049 fSpacePointData[counter].fX = fCurrentRow;
1050 fSpacePointData[counter].fY = fpad;
1051 fSpacePointData[counter].fZ = ftime;
1052 }
1053
1054 fSpacePointData[counter].fCharge = list[j].fTotalCharge;
1055 fSpacePointData[counter].fPadRow = fCurrentRow;
1056 fSpacePointData[counter].fSigmaY2 = fpad2;
1057 fSpacePointData[counter].fSigmaZ2 = ftime2;
1058
6b15c309 1059 fSpacePointData[counter].fMaxQ = list[j].fQMax;
1060
01f43166 1061 fSpacePointData[counter].fUsed = kFALSE; // only used / set in AliHLTTPCDisplay
1062 fSpacePointData[counter].fTrackN = -1; // only used / set in AliHLTTPCDisplay
1063
1064 Int_t patch=fCurrentPatch;
1065 if(patch==-1) patch=0; //never store negative patch number
1066 fSpacePointData[counter].fID = counter
1067 +((fCurrentSlice&0x7f)<<25)+((patch&0x7)<<22);//Uli
1068
1069#ifdef do_mc
1070 Int_t trackID[3];
1071 GetTrackID((Int_t)rint(fpad),(Int_t)rint(ftime),trackID);
1072
1073 fSpacePointData[counter].fTrackID[0] = trackID[0];
1074 fSpacePointData[counter].fTrackID[1] = trackID[1];
1075 fSpacePointData[counter].fTrackID[2] = trackID[2];
1076
01f43166 1077#endif
1078
1079 fNClusters++;
1080 counter++;
1081 }
1082}