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