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