]> 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
Update master to aliroot
[u/mrichter/AliRoot.git] / HLT / TPCLib / AliHLTTPCClusterFinder.cxx
CommitLineData
5dcddb3e 1// $Id$
a38a7850 2
297174de 3//**************************************************************************
4//* This file is property of and copyright by the ALICE HLT Project *
5//* ALICE Experiment at CERN, All rights reserved. *
6//* *
7//* Primary Authors: Anders Vestbo, Constantin Loizides *
8//* Developers: Kenneth Aamodt <kenneth.aamodt@student.uib.no> *
9//* Kalliopi Kanaki *
10//* for The ALICE HLT Project. *
11//* *
12//* Permission to use, copy, modify and distribute this software and its *
13//* documentation strictly for non-commercial purposes is hereby granted *
14//* without fee, provided that the above copyright notice appears in all *
15//* copies and that both the copyright notice and this permission notice *
16//* appear in the supporting documentation. The authors make no claims *
17//* about the suitability of this software for any purpose. It is *
18//* provided "as is" without express or implied warranty. *
19//**************************************************************************
a38a7850 20
b783c3e8 21// @file AliHLTTPCClusterFinder.cxx
22// @author Kenneth Aamodt, Kalliopi Kanaki
23// @date
24// @brief Cluster Finder for the TPC
25// @note
a38a7850 26
27#include "AliHLTTPCDigitReader.h"
a38a7850 28#include "AliHLTTPCLogging.h"
29#include "AliHLTTPCClusterFinder.h"
a38a7850 30#include "AliHLTTPCSpacePointData.h"
46b33a24 31#include "AliHLTTPCPad.h"
8252a538 32#include <sys/time.h>
deba5d85 33#include <algorithm>
6a5788b9 34#include <cmath>
5ee8e602 35#include "AliTPCcalibDB.h"
36#include "AliTPCTransform.h"
a442ba54 37#include "AliTPCParam.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),
6b15c309 53 fDeconvTime(kFALSE),
54 fDeconvPad(kFALSE),
46b33a24 55 fStdout(kFALSE),
56 fCalcerr(kTRUE),
5e75f4e0 57 fFillRawClusters(kFALSE),
46b33a24 58 fFirstRow(0),
59 fLastRow(0),
2a083ac4 60 fCurrentRow(0),
61 fCurrentSlice(0),
62 fCurrentPatch(0),
63 fMatch(1),
64 fThreshold(10),
2a083ac4 65 fNClusters(0),
66 fMaxNClusters(0),
67 fXYErr(0.2),
68 fZErr(0.3),
01f43166 69 fOccupancyLimit(1.0),
64defa03 70 fUnsorted(0),
a74855c2 71 fVectorInitialized(kFALSE),
8252a538 72 fClusters(),
deba5d85 73 fClustersMCInfo(),
74 fMCDigits(),
5e75f4e0 75 fRawClusters(),
a74855c2 76 fNumberOfPadsInRow(NULL),
77 fNumberOfRows(0),
2fdb1ae7 78 fRowOfFirstCandidate(0),
6b15c309 79 fDoPadSelection(kFALSE),
80 fFirstTimeBin(0),
81 fLastTimeBin(AliHLTTPCTransform::GetNTimeBins()),
82 fTotalChargeOfPreviousClusterCandidate(0),
deba5d85 83 fChargeOfCandidatesFalling(kFALSE),
84 f32BitFormat(kFALSE),
85 fDoMC(kFALSE),
5ee8e602 86 fClusterMCVector(),
6a5788b9 87 fOfflineTransform(NULL),
a442ba54 88 fOfflineTPCParam( NULL ),
89 fOfflineTPCRecoParam(*AliTPCRecoParam::GetHLTParam()),
db76ab35 90 fTimeMeanDiff(2),
91 fReleaseMemory(0)
a38a7850 92{
8252a538 93 //constructor
a442ba54 94
95 //uptate the transform class
96
5ee8e602 97 fOfflineTransform = AliTPCcalibDB::Instance()->GetTransform();
98 if(!fOfflineTransform){
a442ba54 99 HLTError("AliHLTTPCClusterFinder()::UpdateCAlibDB:: Offline transform not in AliTPCcalibDB.");
5ee8e602 100 }
f14729cb 101 else{
f14729cb 102 fOfflineTransform->SetCurrentRecoParam(&fOfflineTPCRecoParam);
103 }
a442ba54 104
105 fOfflineTPCParam = AliTPCcalibDB::Instance()->GetParameters();
106 if( !fOfflineTPCParam ){
107 HLTError("AliHLTTPCClusterFinder()::UpdateCAlibDB:: Offline TPC parameters not in AliTPCcalibDB.");
108 } else {
109 fOfflineTPCParam->Update();
110 fOfflineTPCParam->ReadGeoMatrices();
111 }
112
46b33a24 113}
a38a7850 114
a912b63b 115AliHLTTPCClusterFinder::~AliHLTTPCClusterFinder(){
116 // see header file for class documentation
117
a38a7850 118 //destructor
8252a538 119 if(fVectorInitialized){
120 DeInitializePadArray();
121 }
122 if(fNumberOfPadsInRow){
123 delete [] fNumberOfPadsInRow;
124 fNumberOfPadsInRow=NULL;
125 }
a38a7850 126}
a912b63b 127
128void AliHLTTPCClusterFinder::InitSlice(Int_t slice,Int_t patch,Int_t nmaxpoints){
129 // see header file for class documentation
130
a38a7850 131 //init slice
132 fNClusters = 0;
133 fMaxNClusters = nmaxpoints;
134 fCurrentSlice = slice;
135 fCurrentPatch = patch;
a912b63b 136 fFirstRow=AliHLTTPCTransform::GetFirstRow(patch);
137 fLastRow=AliHLTTPCTransform::GetLastRow(patch);
37c9ee0f 138
139 fClusters.clear();
140 fClustersMCInfo.clear();
141 fMCDigits.clear();
142 fClusterMCVector.clear();
5e75f4e0 143 fRawClusters.clear();
a38a7850 144}
145
8252a538 146void AliHLTTPCClusterFinder::InitializePadArray(){
147 // see header file for class documentation
148
149 if(fCurrentPatch>5||fCurrentPatch<0){
150 HLTFatal("Patch is not set");
151 return;
152 }
153
154 HLTDebug("Patch number=%d",fCurrentPatch);
155
156 fFirstRow = AliHLTTPCTransform::GetFirstRow(fCurrentPatch);
157 fLastRow = AliHLTTPCTransform::GetLastRow(fCurrentPatch);
158
159 fNumberOfRows=fLastRow-fFirstRow+1;
160 fNumberOfPadsInRow= new UInt_t[fNumberOfRows];
161
162 memset( fNumberOfPadsInRow, 0, sizeof(Int_t)*(fNumberOfRows));
163
db76ab35 164 fRowPadVector.clear();
165
a74855c2 166 for(UInt_t i=0;i<fNumberOfRows;i++){
8252a538 167 fNumberOfPadsInRow[i]=AliHLTTPCTransform::GetNPads(i+fFirstRow);
168 AliHLTTPCPadVector tmpRow;
288f1d48 169 for(UInt_t j=0;j<=fNumberOfPadsInRow[i];j++){
8252a538 170 AliHLTTPCPad *tmpPad = new AliHLTTPCPad(2);
171 tmpPad->SetID(i,j);
172 tmpRow.push_back(tmpPad);
173 }
174 fRowPadVector.push_back(tmpRow);
175 }
176 fVectorInitialized=kTRUE;
177}
178
a912b63b 179Int_t AliHLTTPCClusterFinder::DeInitializePadArray(){
8252a538 180 // see header file for class documentation
a912b63b 181
db76ab35 182 if( fVectorInitialized ){
183 for(UInt_t i=0;i<fNumberOfRows;i++){
184 for(UInt_t j=0;j<=fNumberOfPadsInRow[i];j++){
185 delete fRowPadVector[i][j];
186 fRowPadVector[i][j]=NULL;
187 }
188 fRowPadVector[i].clear();
8252a538 189 }
db76ab35 190 fRowPadVector.clear();
191 delete[] fNumberOfPadsInRow;
192 fNumberOfPadsInRow = 0;
8252a538 193 }
db76ab35 194 fVectorInitialized=kFALSE;
8252a538 195 return 1;
db76ab35 196}
8252a538 197
a38a7850 198
a912b63b 199void AliHLTTPCClusterFinder::SetOutputArray(AliHLTTPCSpacePointData *pt){
200 // see header file for class documentation
a38a7850 201 //set pointer to output
202 fSpacePointData = pt;
203}
204
a912b63b 205
206void AliHLTTPCClusterFinder::ReadDataUnsorted(void* ptr,unsigned long size){
207 // see header file for class documentation
a38a7850 208 //set input pointer
209 fPtr = (UChar_t*)ptr;
210 fSize = size;
a38a7850 211
a912b63b 212 if(!fVectorInitialized){
213 InitializePadArray();
6b15c309 214 }
a912b63b 215
216 if (fDigitReader->InitBlock(fPtr,fSize,fFirstRow,fLastRow,fCurrentPatch,fCurrentSlice)<0) {
217 HLTError("failed setting up digit reader (InitBlock)");
218 return;
84645eb0 219 }
a912b63b 220
221 while(fDigitReader->NextChannel()){
222 UInt_t row=fDigitReader->GetRow();
223 UInt_t pad=fDigitReader->GetPad();
a38a7850 224
a912b63b 225 if(row>=fRowPadVector.size()){
226 HLTError("Row number is to large: %d, max is %d",row,fRowPadVector.size()-1);
227 continue;
a38a7850 228 }
a912b63b 229 if(pad>=fRowPadVector[row].size()){
230 HLTError("Pad number is to large: %d, max is %d",pad,fRowPadVector[row].size());
231 continue;
a38a7850 232 }
233
a912b63b 234 while(fDigitReader->NextBunch()){
235 if(fDigitReader->GetBunchSize()>1){//to remove single timebin values, this will have to change at some point
236 UInt_t time = fDigitReader->GetTime();
237 if((Int_t)time>=fFirstTimeBin && (Int_t)time+fDigitReader->GetBunchSize()<=fLastTimeBin){
deba5d85 238 // Kenneth: 20-04-09. The following if have been added because of inconsistency in the 40 bit decoder and the 32 bit decoder.
239 // GetSignals() in the 40 bit decoder returns an array of UInt_t while the 32 bit one returns UShort_t
240 // The same is true for the function ReadDataUnsortedDeconvoluteTime() below.
241 // In addition the signals are organized in the opposite direction
242 if(f32BitFormat){
243 const UShort_t *bunchData= fDigitReader->GetSignalsShort();
244 AliHLTTPCClusters candidate;
245 for(Int_t i=fDigitReader->GetBunchSize()-1;i>=0;i--){
246 candidate.fTotalCharge+=bunchData[i];
247 candidate.fTime += time*bunchData[i];
248 candidate.fTime2 += time*time*bunchData[i];
249 if(bunchData[i]>candidate.fQMax){
250 candidate.fQMax=bunchData[i];
251 }
252 time++;
253 }
254 if(candidate.fTotalCharge>0){
255 candidate.fMean=candidate.fTime/candidate.fTotalCharge;
256 candidate.fPad=candidate.fTotalCharge*pad;
257 candidate.fPad2=candidate.fPad*pad;
258 candidate.fLastMergedPad=pad;
259 candidate.fRowNumber=row+fDigitReader->GetRowOffset();
260 }
261 if(fRowPadVector[row][pad] != NULL){
262 fRowPadVector[row][pad]->AddClusterCandidate(candidate);
a912b63b 263 }
46b33a24 264 }
deba5d85 265 else{
266 const UInt_t *bunchData= fDigitReader->GetSignals();
267 AliHLTTPCClusters candidate;
025ba4dc 268 const AliHLTTPCDigitData* digits = NULL;
269 if(fDoMC && (digits = fDigitReader->GetBunchDigits())!=NULL){
deba5d85 270 for(Int_t i=0;i<fDigitReader->GetBunchSize();i++){
271 candidate.fTotalCharge+=bunchData[i];
272 candidate.fTime += time*bunchData[i];
273 candidate.fTime2 += time*time*bunchData[i];
274 if(bunchData[i]>candidate.fQMax){
275 candidate.fQMax=bunchData[i];
276 }
277 fMCDigits.push_back(digits[i]);
278 time++;
279 }
280 }
281 else{
282 for(Int_t i=0;i<fDigitReader->GetBunchSize();i++){
283 candidate.fTotalCharge+=bunchData[i];
284 candidate.fTime += time*bunchData[i];
285 candidate.fTime2 += time*time*bunchData[i];
286 if(bunchData[i]>candidate.fQMax){
287 candidate.fQMax=bunchData[i];
288 }
289 time++;
290 }
291 }
292 if(candidate.fTotalCharge>0){
293 candidate.fMean=candidate.fTime/candidate.fTotalCharge;
294 candidate.fPad=candidate.fTotalCharge*pad;
295 candidate.fPad2=candidate.fPad*pad;
296 candidate.fLastMergedPad=pad;
297 candidate.fRowNumber=row+fDigitReader->GetRowOffset();
298 }
299 if(fRowPadVector[row][pad] != NULL){
300 fRowPadVector[row][pad]->AddClusterCandidate(candidate);
301 if(fDoMC){
302 fRowPadVector[row][pad]->AddCandidateDigits(fMCDigits);
303 fMCDigits.clear();
304 }
305 }
46b33a24 306 }
307 }
308 }
a912b63b 309 }
310 }
311}
a38a7850 312
a912b63b 313void AliHLTTPCClusterFinder::ReadDataUnsortedDeconvoluteTime(void* ptr,unsigned long size){
314 // see header file for class documentation
738c049f 315
a912b63b 316 //set input pointer
317 fPtr = (UChar_t*)ptr;
318 fSize = size;
738c049f 319
a912b63b 320 if(!fVectorInitialized){
321 InitializePadArray();
322 }
738c049f 323
a912b63b 324 if (fDigitReader->InitBlock(fPtr,fSize,fFirstRow,fLastRow,fCurrentPatch,fCurrentSlice)<0) {
325 HLTError("failed setting up digit reader (InitBlock)");
326 return;
327 }
328
329 while(fDigitReader->NextChannel()){
330 UInt_t row=fDigitReader->GetRow();
331 UInt_t pad=fDigitReader->GetPad();
46b33a24 332
a912b63b 333 while(fDigitReader->NextBunch()){
334 if(fDigitReader->GetBunchSize()>1){//to remove single timebin values, this will have to change at some point
335 UInt_t time = fDigitReader->GetTime();
336 if((Int_t)time>=fFirstTimeBin && (Int_t)time+fDigitReader->GetBunchSize()<=fLastTimeBin){
337 Int_t indexInBunchData=0;
338 Bool_t moreDataInBunch=kFALSE;
339 UInt_t prevSignal=0;
340 Bool_t signalFalling=kFALSE;
deba5d85 341
342 // Kenneth: 20-04-09. The following if have been added because of inconsistency in the 40 bit decoder and the 32 bit decoder.
343 // GetSignals() in the 40 bit decoder returns an array of UInt_t while the 32 bit one returns UShort_t
344 // The same is true for the function ReadDataUnsorted() above.
345 // In addition the signals are organized in the opposite direction
346 if(f32BitFormat){
256769f4 347 indexInBunchData = fDigitReader->GetBunchSize()-1;
deba5d85 348 const UShort_t *bunchData= fDigitReader->GetSignalsShort();
256769f4 349
deba5d85 350 do{
351 AliHLTTPCClusters candidate;
352 //for(Int_t i=indexInBunchData;i<fDigitReader->GetBunchSize();i++){
353 for(Int_t i=indexInBunchData;i>=0;i--){
256769f4 354
deba5d85 355 // Checks if one need to deconvolute the signals
356 if(bunchData[i]>prevSignal && signalFalling==kTRUE){
357 if(i<fDigitReader->GetBunchSize()-1){ // means there are more than one signal left in the bunch
358 moreDataInBunch=kTRUE;
359 prevSignal=0;
360 }
361 break;
362 }
363
364 // Checks if the signal is 0, then quit processing the data.
365 if(bunchData[i]==0 && i<fDigitReader->GetBunchSize()-1){//means we have 0 data fom the rcu, might happen depending on the configuration settings
a912b63b 366 moreDataInBunch=kTRUE;
367 prevSignal=0;
deba5d85 368 break;
369 }
370
371 if(prevSignal>bunchData[i]){//means the peak of the signal has been reached and deconvolution will happen if the signal rise again.
372 signalFalling=kTRUE;
373 }
374 candidate.fTotalCharge+=bunchData[i];
375 candidate.fTime += time*bunchData[i];
376 candidate.fTime2 += time*time*bunchData[i];
377 if(bunchData[i]>candidate.fQMax){
378 candidate.fQMax=bunchData[i];
a912b63b 379 }
deba5d85 380 prevSignal=bunchData[i];
381 time++;
382 indexInBunchData--;
a912b63b 383 }
deba5d85 384 if(candidate.fTotalCharge>0){
385 candidate.fMean=candidate.fTime/candidate.fTotalCharge;
386 candidate.fPad=candidate.fTotalCharge*pad;
387 candidate.fPad2=candidate.fPad*pad;
388 candidate.fLastMergedPad=pad;
389 candidate.fRowNumber=row+fDigitReader->GetRowOffset();
a912b63b 390 }
0f98f30d 391 fRowPadVector[row][pad]->AddClusterCandidate(candidate);
deba5d85 392 if(indexInBunchData<fDigitReader->GetBunchSize()-1){
393 moreDataInBunch=kFALSE;
a912b63b 394 }
deba5d85 395 }while(moreDataInBunch);
396 }
397 else{
398 const UInt_t *bunchData= fDigitReader->GetSignals();
399 do{
400 AliHLTTPCClusters candidate;
0f98f30d 401 const AliHLTTPCDigitData* digits = fDigitReader->GetBunchDigits();
402 if(fDoMC) fMCDigits.clear();
403
deba5d85 404 for(Int_t i=indexInBunchData;i<fDigitReader->GetBunchSize();i++){
405 // Checks if one need to deconvolute the signals
406 if(bunchData[i]>prevSignal && signalFalling==kTRUE){
407 if(i<fDigitReader->GetBunchSize()-1){ // means there are more than one signal left in the bunch
408 moreDataInBunch=kTRUE;
409 prevSignal=0;
410 }
411 break;
412 }
413
414 // Checks if the signal is 0, then quit processing the data.
415 if(bunchData[i]==0 && i<fDigitReader->GetBunchSize()-1){//means we have 0 data fom the rcu, might happen depending on the configuration settings
416 moreDataInBunch=kTRUE;
417 prevSignal=0;
418 break;
419 }
420
421 if(prevSignal>bunchData[i]){//means the peak of the signal has been reached and deconvolution will happen if the signal rise again.
422 signalFalling=kTRUE;
423 }
424 candidate.fTotalCharge+=bunchData[i];
425 candidate.fTime += time*bunchData[i];
426 candidate.fTime2 += time*time*bunchData[i];
427 if(bunchData[i]>candidate.fQMax){
428 candidate.fQMax=bunchData[i];
429 }
0f98f30d 430 if( fDoMC ) fMCDigits.push_back(digits[i]);
deba5d85 431
432 prevSignal=bunchData[i];
433 time++;
434 indexInBunchData++;
a912b63b 435 }
deba5d85 436 if(candidate.fTotalCharge>0){
437 candidate.fMean=candidate.fTime/candidate.fTotalCharge;
438 candidate.fPad=candidate.fTotalCharge*pad;
439 candidate.fPad2=candidate.fPad*pad;
440 candidate.fLastMergedPad=pad;
441 candidate.fRowNumber=row+fDigitReader->GetRowOffset();
442 }
443 fRowPadVector[row][pad]->AddClusterCandidate(candidate);
0f98f30d 444 if(fDoMC){
445 fRowPadVector[row][pad]->AddCandidateDigits(fMCDigits);
446 fMCDigits.clear();
447 }
deba5d85 448 if(indexInBunchData<fDigitReader->GetBunchSize()-1){
449 moreDataInBunch=kFALSE;
450 }
451 }while(moreDataInBunch);
452 }
a912b63b 453 }
454 }
455 }
456 }
457}
458
459Bool_t AliHLTTPCClusterFinder::ComparePads(AliHLTTPCPad *nextPad,AliHLTTPCClusters* cluster,Int_t nextPadToRead){
460 // see header file for class documentation
461
462 //Checking if we have a match on the next pad
463 for(UInt_t candidateNumber=0;candidateNumber<nextPad->fClusterCandidates.size();candidateNumber++){
2d6b17ff 464 if(nextPad->fUsedClusterCandidates[candidateNumber] == 1){
465 continue;
466 }
a912b63b 467 AliHLTTPCClusters *candidate =&nextPad->fClusterCandidates[candidateNumber];
6a5788b9 468 // if(cluster->fMean-candidate->fMean==1 || candidate->fMean-cluster->fMean==1 || cluster->fMean-candidate->fMean==0){
469
470 if( abs((Int_t)(cluster->fMean - candidate->fMean)) <= fTimeMeanDiff ){
a912b63b 471 if(fDeconvPad){
472 if(candidate->fTotalCharge<fTotalChargeOfPreviousClusterCandidate){//peak is reached
473 fChargeOfCandidatesFalling=kTRUE;
474 }
475 if(candidate->fTotalCharge>fTotalChargeOfPreviousClusterCandidate && fChargeOfCandidatesFalling==kTRUE){//we have deconvolution
476 return kFALSE;
477 }
478 }
479 cluster->fMean=candidate->fMean;
480 cluster->fTotalCharge+=candidate->fTotalCharge;
481 cluster->fTime += candidate->fTime;
482 cluster->fTime2 += candidate->fTime2;
483 cluster->fPad+=candidate->fPad;
484 cluster->fPad2+=candidate->fPad2;
485 cluster->fLastMergedPad=candidate->fPad;
486 if(candidate->fQMax>cluster->fQMax){
487 cluster->fQMax=candidate->fQMax;
488 }
deba5d85 489 if(fDoMC){
490 FillMCClusterVector(nextPad->GetCandidateDigits(candidateNumber));
491 }
492
a912b63b 493 if(fDoPadSelection){
494 UInt_t rowNo = nextPad->GetRowNumber();
495 UInt_t padNo = nextPad->GetPadNumber();
496 if(padNo-1>0){
497 fRowPadVector[rowNo][padNo-2]->fSelectedPad=kTRUE;
498 fRowPadVector[rowNo][padNo-2]->fHWAddress=(AliHLTUInt16_t)fDigitReader->GetAltroBlockHWaddr(rowNo,padNo-2);
499 }
500 fRowPadVector[rowNo][padNo-1]->fSelectedPad=kTRUE;// quick solution to set the first pad to selected
501 fRowPadVector[rowNo][padNo-1]->fHWAddress=(AliHLTUInt16_t)fDigitReader->GetAltroBlockHWaddr(rowNo,padNo-1);
502 fRowPadVector[rowNo][padNo]->fSelectedPad=kTRUE;
503 fRowPadVector[rowNo][padNo]->fHWAddress=(AliHLTUInt16_t)fDigitReader->GetAltroBlockHWaddr(rowNo,padNo);
504 }
505
506 //setting the matched pad to used
507 nextPad->fUsedClusterCandidates[candidateNumber]=1;
508 nextPadToRead++;
509 if(nextPadToRead<(Int_t)fNumberOfPadsInRow[fRowOfFirstCandidate]){
510 nextPad=fRowPadVector[fRowOfFirstCandidate][nextPadToRead];
511 ComparePads(nextPad,cluster,nextPadToRead);
512 }
513 else{
514 return kFALSE;
515 }
516 }
a912b63b 517 }
518 return kFALSE;
519}
520
521Int_t AliHLTTPCClusterFinder::FillHWAddressList(AliHLTUInt16_t *hwaddlist, Int_t maxHWadd){
522 // see header file for class documentation
523
524 Int_t counter=0;
525 for(UInt_t row=0;row<fNumberOfRows;row++){
526 for(UInt_t pad=0;pad<fNumberOfPadsInRow[row]-1;pad++){
527 if(fRowPadVector[row][pad]->fSelectedPad){
528 if(counter<maxHWadd){
529 hwaddlist[counter]=(AliHLTUInt16_t)fRowPadVector[row][pad]->fHWAddress;
530 counter++;
531 }
532 else{
533 HLTWarning("To many hardwareaddresses, skip adding");
534 }
535
536 }
537 }
538 }
539 return counter;
540}
deba5d85 541
542
d29671e2 543Int_t AliHLTTPCClusterFinder::FillOutputMCInfo(AliHLTTPCClusterMCLabel * outputMCInfo, Int_t maxNumberOfClusterMCInfo){
deba5d85 544 // see header file for class documentation
545
546 Int_t counter=0;
547
548 for(UInt_t mc=0;mc<fClustersMCInfo.size();mc++){
549 if(counter<maxNumberOfClusterMCInfo){
550 outputMCInfo[counter] = fClustersMCInfo[mc];
551 counter++;
552 }
553 else{
554 HLTWarning("To much MCInfo has been added (no more space), skip adding");
555 }
556 }
557 return counter;
558}
a912b63b 559
5e75f4e0 560Int_t AliHLTTPCClusterFinder::FillOutputRaw(AliHLTTPCRawCluster* rawClusters, unsigned sizeInByte) const
561{
562 // fill the raw clusters
563 if (fRawClusters.size()*sizeof(AliHLTTPCRawCluster)>sizeInByte) {
564 HLTError("not enough space to write raw clusters");
565 return 0;
566 }
567 memcpy(rawClusters, &fRawClusters[0], fRawClusters.size()*sizeof(AliHLTTPCRawCluster));
568 return fRawClusters.size();
569}
570
a912b63b 571void AliHLTTPCClusterFinder::FindClusters(){
572 // see header file for function documentation
573
574 AliHLTTPCClusters* tmpCandidate=NULL;
575 for(UInt_t row=0;row<fNumberOfRows;row++){
576 fRowOfFirstCandidate=row;
577 for(UInt_t pad=0;pad<fNumberOfPadsInRow[row];pad++){
578 AliHLTTPCPad *tmpPad=fRowPadVector[row][pad];
579 for(size_t candidate=0;candidate<tmpPad->fClusterCandidates.size();candidate++){
580 if(tmpPad->fUsedClusterCandidates[candidate]){
581 continue;
582 }
583 tmpCandidate=&tmpPad->fClusterCandidates[candidate];
584 UInt_t tmpTotalCharge=tmpCandidate->fTotalCharge;
585
deba5d85 586 if(fDoMC){
587 fClusterMCVector.clear();
588 FillMCClusterVector(tmpPad->GetCandidateDigits(candidate));
589 }
590
a912b63b 591 ComparePads(fRowPadVector[row][pad+1],tmpCandidate,pad+1);
592 if(tmpCandidate->fTotalCharge>tmpTotalCharge){
593 //we have a cluster
594 fClusters.push_back(*tmpCandidate);
deba5d85 595 if(fDoMC){
596 //sort the vector (large->small) according to weight and remove elements above 2 (keep 0 1 and 2)
d29671e2 597 sort(fClusterMCVector.begin(),fClusterMCVector.end(), CompareWeights );
598 AliHLTTPCClusterMCLabel tmpClusterMCInfo;
deba5d85 599
d29671e2 600 AliHLTTPCClusterMCWeight zeroMC;
deba5d85 601 zeroMC.fMCID=-1;
602 zeroMC.fWeight=0;
603
604 if(fClusterMCVector.size()>0){
605 tmpClusterMCInfo.fClusterID[0]=fClusterMCVector.at(0);
606 }
607 else{
608 tmpClusterMCInfo.fClusterID[0]=zeroMC;
609 }
610
611 if(fClusterMCVector.size()>1){
612 tmpClusterMCInfo.fClusterID[1]=fClusterMCVector.at(1);
613 }
614 else{
615 tmpClusterMCInfo.fClusterID[1]=zeroMC;
616 }
617
618 if(fClusterMCVector.size()>2){
619 tmpClusterMCInfo.fClusterID[2]=fClusterMCVector.at(2);
620 }
621 else{
622 tmpClusterMCInfo.fClusterID[2]=zeroMC;
623 }
624
625 fClustersMCInfo.push_back(tmpClusterMCInfo);
626 }
627
a912b63b 628 }
629 }
630 tmpPad->ClearCandidates();
631 }
632 fRowPadVector[row][fNumberOfPadsInRow[row]]->ClearCandidates();
633 }
634
635 HLTInfo("Found %d clusters.",fClusters.size());
636
637 //TODO: Change so it stores AliHLTTPCSpacePointData directly, instead of this copying
638
639 AliClusterData * clusterlist = new AliClusterData[fClusters.size()]; //Clusterlist
640 for(unsigned int i=0;i<fClusters.size();i++){
641 clusterlist[i].fTotalCharge = fClusters[i].fTotalCharge;
642 clusterlist[i].fPad = fClusters[i].fPad;
643 clusterlist[i].fPad2 = fClusters[i].fPad2;
644 clusterlist[i].fTime = fClusters[i].fTime;
645 clusterlist[i].fTime2 = fClusters[i].fTime2;
646 clusterlist[i].fMean = fClusters[i].fMean;
647 clusterlist[i].fFlags = fClusters[i].fFlags;
648 clusterlist[i].fChargeFalling = fClusters[i].fChargeFalling;
649 clusterlist[i].fLastCharge = fClusters[i].fLastCharge;
650 clusterlist[i].fLastMergedPad = fClusters[i].fLastMergedPad;
651 clusterlist[i].fRow = fClusters[i].fRowNumber;
652 clusterlist[i].fQMax = fClusters[i].fQMax;
653 }
654
655 WriteClusters(fClusters.size(),clusterlist);
656 delete [] clusterlist;
657 fClusters.clear();
db76ab35 658 if( fReleaseMemory ) DeInitializePadArray();// call this when the -releaseMemory flag is set
a912b63b 659}
660
661
4f43db26 662Bool_t AliHLTTPCClusterFinder::UpdateCalibDB(){
663
664 //update the db
665 AliTPCcalibDB::Instance()->Update();
666
a442ba54 667 Bool_t ret = 1;
668
4f43db26 669 //uptate the transform class
a442ba54 670
671 fOfflineTransform = AliTPCcalibDB::Instance()->GetTransform();
672 if(!fOfflineTransform){
673 HLTError("AliHLTTPCClusterFinder()::UpdateCAlibDB:: Offline transform not in AliTPCcalibDB.");
674 ret = 0;
f4f7aad5 675 }
a442ba54 676 else{
677 fOfflineTransform->SetCurrentRecoParam(&fOfflineTPCRecoParam);
678 }
679
680 fOfflineTPCParam = AliTPCcalibDB::Instance()->GetParameters();
681 if( !fOfflineTPCParam ){
682 HLTError("AliHLTTPCClusterFinder()::UpdateCAlibDB:: Offline TPC parameters not in AliTPCcalibDB.");
683 ret = 0;
684 } else {
685 fOfflineTPCParam->Update();
686 fOfflineTPCParam->ReadGeoMatrices();
687 }
688
689 return ret;
4f43db26 690}
691
a912b63b 692//---------------------------------- Under this line the old sorted clusterfinder functions can be found --------------------------------
693
694
695void AliHLTTPCClusterFinder::PrintClusters(){
696 // see header file for class documentation
697
698 for(size_t i=0;i<fClusters.size();i++){
699 HLTInfo("Cluster number: %d",i);
700 HLTInfo("Row: %d \t Pad: %d",fClusters[i].fRowNumber,fClusters[i].fPad/fClusters[i].fTotalCharge);
701 HLTInfo("Total Charge: %d",fClusters[i].fTotalCharge);
702 HLTInfo("fPad: %d",fClusters[i].fPad);
703 HLTInfo("PadError: %d",fClusters[i].fPad2);
704 HLTInfo("TimeMean: %d",fClusters[i].fTime/fClusters[i].fTotalCharge);
705 HLTInfo("TimeError: %d",fClusters[i].fTime2);
706 HLTInfo("EndOfCluster:");
707 }
708}
709
0f98f30d 710void AliHLTTPCClusterFinder::FillMCClusterVector(vector<AliHLTTPCDigitData> *digitData){
b783c3e8 711 // see header file for class documentation
0f98f30d 712 if( !digitData ) return;
713 for(UInt_t d=0;d<digitData->size();d++){
714 Int_t nIDsInDigit = (digitData->at(d).fTrackID[0]>=0) + (digitData->at(d).fTrackID[1]>=0) + (digitData->at(d).fTrackID[2]>=0);
deba5d85 715 for(Int_t id=0; id<3; id++){
0f98f30d 716 if(digitData->at(d).fTrackID[id]>=0){
deba5d85 717 Bool_t matchFound = kFALSE;
d29671e2 718 AliHLTTPCClusterMCWeight mc;
0f98f30d 719 mc.fMCID = digitData->at(d).fTrackID[id];
720 mc.fWeight = ((Float_t)digitData->at(d).fCharge)/nIDsInDigit;
deba5d85 721 for(UInt_t i=0;i<fClusterMCVector.size();i++){
722 if(mc.fMCID == fClusterMCVector.at(i).fMCID){
723 fClusterMCVector.at(i).fWeight += mc.fWeight;
724 matchFound = kTRUE;
725 }
726 }
727 if(matchFound == kFALSE){
728 fClusterMCVector.push_back(mc);
729 }
730 }
731 }
732 }
733}
734
735
a912b63b 736void AliHLTTPCClusterFinder::Read(void* ptr,unsigned long size){
737 //set input pointer
738 fPtr = (UChar_t*)ptr;
739 fSize = size;
740}
741
742void AliHLTTPCClusterFinder::ProcessDigits(){
743 // see header file for class documentation
744
745 int iResult=0;
746 bool readValue = true;
747 Int_t newRow = 0;
748 Int_t rowOffset = 0;
749 UShort_t time=0,newTime=0;
750 UInt_t pad=0,newPad=0;
751 AliHLTTPCSignal_t charge=0;
752
753 fNClusters = 0;
754
755 // initialize block for reading packed data
756 iResult=fDigitReader->InitBlock(fPtr,fSize,fFirstRow,fLastRow,fCurrentPatch,fCurrentSlice);
757 if (iResult<0) return;
758
759 readValue = fDigitReader->Next();
760
761 // Matthias 08.11.2006 the following return would cause termination without writing the
762 // ClusterData and thus would block the component. I just want to have the commented line
763 // here for information
764 //if (!readValue)return;
765
766 pad = fDigitReader->GetPad();
767 time = fDigitReader->GetTime();
768 fCurrentRow = fDigitReader->GetRow();
769
770 if ( fCurrentPatch >= 2 ) // Outer sector, patches 2, 3, 4, 5
771 rowOffset = AliHLTTPCTransform::GetFirstRow( 2 );
772
773 fCurrentRow += rowOffset;
774
775 UInt_t lastpad = 123456789;
776 const UInt_t kPadArraySize=5000;
777 const UInt_t kClusterListSize=10000;
778 AliClusterData *pad1[kPadArraySize]; //2 lists for internal memory=2pads
779 AliClusterData *pad2[kPadArraySize]; //2 lists for internal memory=2pads
780 AliClusterData clusterlist[kClusterListSize]; //Clusterlist
781
782 AliClusterData **currentPt; //List of pointers to the current pad
783 AliClusterData **previousPt; //List of pointers to the previous pad
784 currentPt = pad2;
785 previousPt = pad1;
786 UInt_t nprevious=0,ncurrent=0,ntotal=0;
787
788 /* quick implementation of baseline calculation and zero suppression
789 open a pad object for each pad and delete it after processing.
790 later a list of pad objects with base line history can be used
791 The whole thing only works if we really get unprocessed raw data, if
792 the data is already zero suppressed, there might be gaps in the time
793 bins.
794 */
795 Int_t gatingGridOffset=50;
796 if(fFirstTimeBin>0){
797 gatingGridOffset=fFirstTimeBin;
798 }
799 AliHLTTPCPad baseline(gatingGridOffset, AliHLTTPCTransform::GetNTimeBins());
800 // just to make later conversion to a list of objects easier
801 AliHLTTPCPad* pCurrentPad=NULL;
802 /*
803 if (fSignalThreshold>=0) {
804 pCurrentPad=&baseline;
805 baseline.SetThreshold(fSignalThreshold);
806 }
807 */
808 while ( readValue!=0 && iResult>=0){ // Reads through all digits in block
809 iResult=0;
810
811 if(pad != lastpad){
812 //This is a new pad
813
814 //Switch the lists:
815 if(currentPt == pad2){
816 currentPt = pad1;
817 previousPt = pad2;
818 }
819 else {
820 currentPt = pad2;
821 previousPt = pad1;
822 }
823 nprevious = ncurrent;
824 ncurrent = 0;
825 if(pad != lastpad+1){
826 //this happens if there is a pad with no signal.
827 nprevious = ncurrent = 0;
828 }
829 lastpad = pad;
830 }
831
832 Bool_t newcluster = kTRUE;
833 UInt_t seqcharge=0,seqaverage=0,seqerror=0;
834 AliHLTTPCSignal_t lastcharge=0;
835 UInt_t bLastWasFalling=0;
836 Int_t newbin=-1;
837
838
839 if(fDeconvTime){
840 redo: //This is a goto.
841
842 if(newbin > -1){
843 //bin = newbin;
844 newbin = -1;
845 }
846
847 lastcharge=0;
848 bLastWasFalling = 0;
849 }
850
851 while(iResult>=0){ //Loop over time bins of current pad
852 // read all the values for one pad at once to calculate the base line
853 if (pCurrentPad) {
854 if (!pCurrentPad->IsStarted()) {
855 //HLTDebug("reading data for pad %d, padrow %d", fDigitReader->GetPad(), fDigitReader->GetRow()+rowOffset);
856 pCurrentPad->SetID(fDigitReader->GetRow()+rowOffset,fDigitReader->GetPad());
857 if ((pCurrentPad->StartEvent())>=0) {
858 do {
859 if ((fDigitReader->GetRow()+rowOffset)!=pCurrentPad->GetRowNumber()) break;
860 if (fDigitReader->GetPad()!=pCurrentPad->GetPadNumber()) break;
861 pCurrentPad->SetRawData(fDigitReader->GetTime(), fDigitReader->GetSignal());
862 //HLTDebug("set raw data to pad: bin %d charge %d", fDigitReader->GetTime(), fDigitReader->GetSignal());
863 } while ((readValue = fDigitReader->Next())!=0);
864 }
865 pCurrentPad->CalculateBaseLine(AliHLTTPCTransform::GetNTimeBins()/2);
866 if (pCurrentPad->Next(kTRUE/*do zero suppression*/)==0) {
867 //HLTDebug("no data available after zero suppression");
868 pCurrentPad->StopEvent();
869 pCurrentPad->ResetHistory();
870 break;
871 }
872 time=pCurrentPad->GetCurrentPosition();
873 if (time>pCurrentPad->GetSize()) {
874 HLTError("invalid time bin for pad");
875 break;
876 }
877 }
878 }
879 if (pCurrentPad) {
880 Float_t occupancy=pCurrentPad->GetOccupancy();
881 //HLTDebug("pad %d occupancy level: %f", pCurrentPad->GetPadNumber(), occupancy);
882 if ( occupancy < fOccupancyLimit ) {
883 charge = pCurrentPad->GetCorrectedData();
884 } else {
885 charge = 0;
886 //HLTDebug("ignoring pad %d with occupancy level %f", pCurrentPad->GetPadNumber(), occupancy);
887 }
888 } else {
889 charge = fDigitReader->GetSignal();
890 }
891 //HLTDebug("get next charge value: position %d charge %d", time, charge);
892
893
894 // CHARGE DEBUG
895 if (fDigitReader->GetRow() == 90){
896///// LOG(AliHLTTPCLog::kFatal,"AliHLTTPCClusterFinder::Row","row90") << "PAD=" << fDigitReader->GetPad() << " TIME=" << fDigitReader->GetTime()
897 // << " SIGNAL=" << fDigitReader->GetSignal() << ENDLOG;
898
899 }
900
901 if(time >= AliHLTTPCTransform::GetNTimeBins()){
902 HLTWarning("Pad %d: Timebin (%d) out of range (%d)", pad, time, AliHLTTPCTransform::GetNTimeBins());
903 iResult=-ERANGE;
904 }
905
906
907 //Get the current ADC-value
a38a7850 908 if(fDeconvTime){
909
910 //Check if the last pixel in the sequence is smaller than this
911 if(charge > lastcharge){
2a083ac4 912 if(bLastWasFalling){
a38a7850 913 newbin = 1;
914 break;
915 }
916 }
2a083ac4 917 else bLastWasFalling = 1; //last pixel was larger than this
a38a7850 918 lastcharge = charge;
919 }
920
921 //Sum the total charge of this sequence
922 seqcharge += charge;
923 seqaverage += time*charge;
924 seqerror += time*time*charge;
925
46b33a24 926 if (pCurrentPad) {
927
928 if((pCurrentPad->Next(kTRUE/*do zero suppression*/))==0) {
929 pCurrentPad->StopEvent();
930 pCurrentPad->ResetHistory();
931 if(readValue) {
932 newPad = fDigitReader->GetPad();
933 newTime = fDigitReader->GetTime();
934 newRow = fDigitReader->GetRow() + rowOffset;
935 }
936 break;
937 }
938
939 newPad=pCurrentPad->GetPadNumber();
940 newTime=pCurrentPad->GetCurrentPosition();
941 newRow=pCurrentPad->GetRowNumber();
942 } else {
a38a7850 943 readValue = fDigitReader->Next();
a38a7850 944 //Check where to stop:
945 if(!readValue) break; //No more value
946
947 newPad = fDigitReader->GetPad();
948 newTime = fDigitReader->GetTime();
949 newRow = fDigitReader->GetRow() + rowOffset;
46b33a24 950 }
a38a7850 951
952 if(newPad != pad)break; //new pad
953 if(newTime != time+1) break; //end of sequence
6b0fdc0e 954 if(iResult<0) break;
a38a7850 955
956 // pad = newpad; is equal
957 time = newTime;
958
959 }//end loop over sequence
960
46b33a24 961 //HLTDebug("ended time bin sequence loop: seqcharge=%d readValue=%d", seqcharge, readValue);
962 //HLTDebug("pad=%d newpad=%d current row=%d newrow=%d", pad, newPad, fCurrentRow, newRow);
963 if (seqcharge<=0) {
964 // with active zero suppression zero values are possible
965 continue;
966 }
a38a7850 967
968 //Calculate mean of sequence:
969 Int_t seqmean=0;
970 if(seqcharge)
971 seqmean = seqaverage/seqcharge;
972 else{
973 LOG(AliHLTTPCLog::kFatal,"AliHLTTPCClusterFinder::ProcessRow","Data")
974 <<"Error in data given to the cluster finder"<<ENDLOG;
975 seqmean = 1;
976 seqcharge = 1;
977 }
978
979 //Calculate mean in pad direction:
980 Int_t padmean = seqcharge*pad;
981 Int_t paderror = pad*padmean;
982
a38a7850 983 //Compare with results on previous pad:
84645eb0 984 for(UInt_t p=0; p<nprevious && p<kPadArraySize && ncurrent<kPadArraySize; p++){
a38a7850 985
986 //dont merge sequences on the same pad twice
987 if(previousPt[p]->fLastMergedPad==pad) continue;
988
989 Int_t difference = seqmean - previousPt[p]->fMean;
990 if(difference < -fMatch) break;
991
992 if(difference <= fMatch){ //There is a match here!!
993 AliClusterData *local = previousPt[p];
994
995 if(fDeconvPad){
996 if(seqcharge > local->fLastCharge){
997 if(local->fChargeFalling){ //The previous pad was falling
998 break; //create a new cluster
999 }
1000 }
1001 else local->fChargeFalling = 1;
1002 local->fLastCharge = seqcharge;
1003 }
1004
1005 //Don't create a new cluster, because we found a match
1006 newcluster = kFALSE;
1007
1008 //Update cluster on current pad with the matching one:
1009 local->fTotalCharge += seqcharge;
1010 local->fPad += padmean;
1011 local->fPad2 += paderror;
1012 local->fTime += seqaverage;
1013 local->fTime2 += seqerror;
1014 local->fMean = seqmean;
1015 local->fFlags++; //means we have more than one pad
1016 local->fLastMergedPad = pad;
1017
1018 currentPt[ncurrent] = local;
1019 ncurrent++;
1020
1021 break;
1022 } //Checking for match at previous pad
1023 } //Loop over results on previous pad.
1024
84645eb0 1025 if(newcluster && ncurrent<kPadArraySize){
a38a7850 1026 //Start a new cluster. Add it to the clusterlist, and update
1027 //the list of pointers to clusters in current pad.
1028 //current pad will be previous pad on next pad.
1029
1030 //Add to the clusterlist:
1031 AliClusterData *tmp = &clusterlist[ntotal];
1032 tmp->fTotalCharge = seqcharge;
1033 tmp->fPad = padmean;
1034 tmp->fPad2 = paderror;
1035 tmp->fTime = seqaverage;
1036 tmp->fTime2 = seqerror;
1037 tmp->fMean = seqmean;
1038 tmp->fFlags = 0; //flags for single pad clusters
1039 tmp->fLastMergedPad = pad;
1040
1041 if(fDeconvPad){
1042 tmp->fChargeFalling = 0;
1043 tmp->fLastCharge = seqcharge;
1044 }
1045
1046 //Update list of pointers to previous pad:
1047 currentPt[ncurrent] = &clusterlist[ntotal];
1048 ntotal++;
1049 ncurrent++;
1050 }
1051
1052 if(fDeconvTime)
1053 if(newbin >= 0) goto redo;
1054
1055 // to prevent endless loop
1056 if(time >= AliHLTTPCTransform::GetNTimeBins()){
46b33a24 1057 HLTWarning("Timebin (%d) out of range (%d)", time, AliHLTTPCTransform::GetNTimeBins());
a38a7850 1058 break;
1059 }
1060
1061
1062 if(!readValue) break; //No more value
84645eb0 1063
1064 if (ntotal>=kClusterListSize || ncurrent>=kPadArraySize) {
1065 HLTWarning("pad array size exceeded ntotal=%d ncurrent=%d, skip rest of the data", ntotal, ncurrent);
1066 break;
1067 }
a38a7850 1068
1069 if(fCurrentRow != newRow){
1070 WriteClusters(ntotal,clusterlist);
1071
1072 lastpad = 123456789;
1073
1074 currentPt = pad2;
1075 previousPt = pad1;
1076 nprevious=0;
1077 ncurrent=0;
1078 ntotal=0;
1079
1080 fCurrentRow = newRow;
1081 }
1082
1083 pad = newPad;
1084 time = newTime;
db16520a 1085
a38a7850 1086 } // END while(readValue)
1087
1088 WriteClusters(ntotal,clusterlist);
1089
46b33a24 1090 HLTInfo("ClusterFinder found %d clusters in slice %d patch %d", fNClusters, fCurrentSlice, fCurrentPatch);
a38a7850 1091
1092} // ENDEND
1093
a912b63b 1094void AliHLTTPCClusterFinder::WriteClusters(Int_t nclusters,AliClusterData *list){
1095 // see header file for class documentation
1096
a38a7850 1097 //write cluster to output pointer
8f8bf0af 1098 Int_t thisrow=-1,thissector=-1;
a38a7850 1099 UInt_t counter = fNClusters;
1100
5e75f4e0 1101 if (fFillRawClusters) {
1102 fRawClusters.resize(nclusters);
1103 }
1104
a38a7850 1105 for(int j=0; j<nclusters; j++)
1106 {
d8b2f74a 1107
1108
1109
deba5d85 1110 if(!list[j].fFlags){
1111 if(fDoMC){
5ee8e602 1112 if(j+(Int_t)fClustersMCInfo.size()-nclusters >=0 && j+fClustersMCInfo.size()-nclusters < fClustersMCInfo.size()){
03efd623 1113 fClustersMCInfo.erase(fClustersMCInfo.begin()+j+fClustersMCInfo.size()-nclusters); // remove the mc info for this cluster since it is not taken into account
1114 }
deba5d85 1115 }
1116 continue; //discard single pad clusters
1117 }
1118 if(list[j].fTotalCharge < fThreshold){
1119 if(fDoMC){
5ee8e602 1120 if(j+(Int_t)fClustersMCInfo.size()-nclusters >=0 && j+fClustersMCInfo.size()-nclusters < fClustersMCInfo.size()){
03efd623 1121 fClustersMCInfo.erase(fClustersMCInfo.begin()+j+fClustersMCInfo.size()-nclusters); // remove the mc info for this cluster since it is not taken into account
1122 }
deba5d85 1123 }
1124 continue; //noise cluster
1125 }
a38a7850 1126 Float_t xyz[3];
1127 Float_t fpad =(Float_t)list[j].fPad / list[j].fTotalCharge;
1128 Float_t fpad2=fXYErr*fXYErr; //fixed given error
1129 Float_t ftime =(Float_t)list[j].fTime / list[j].fTotalCharge;
1130 Float_t ftime2=fZErr*fZErr; //fixed given error
1131
738c049f 1132
aff6e981 1133
a912b63b 1134 if(fUnsorted){
1135 fCurrentRow=list[j].fRow;
8252a538 1136 }
8252a538 1137
a912b63b 1138
1139 if(fCalcerr) { //calc the errors, otherwice take the fixed error
1140 Int_t patch = AliHLTTPCTransform::GetPatch(fCurrentRow);
1141 UInt_t q2=list[j].fTotalCharge*list[j].fTotalCharge;
1142 // Float_t sy2=list[j].fPad2 * list[j].fTotalCharge - list[j].fPad * list[j].fPad;
1143 Float_t sy2=(Float_t)list[j].fPad2 * list[j].fTotalCharge - (Float_t)list[j].fPad * list[j].fPad;
1ee4d015 1144 if(q2 == 0) {
1145 LOG(AliHLTTPCLog::kError,"AliHLTTPCClusterFinder::WriteClusters","Cluster width")
1146 <<"zero charge "<< list[j].fTotalCharge <<" on row "<<fCurrentRow<<" "<<fpad<<" "<<ftime<<ENDLOG;
1147 continue;
1148 }
a912b63b 1149 sy2/=q2;
1150 if(sy2 < 0) {
1151 LOG(AliHLTTPCLog::kError,"AliHLTTPCClusterFinder::WriteClusters","Cluster width")
1152 <<"SigmaY2 negative "<<sy2<<" on row "<<fCurrentRow<<" "<<fpad<<" "<<ftime<<ENDLOG;
1153 continue;
1154 } else {
5e75f4e0 1155 {
a912b63b 1156 fpad2 = (sy2 + 1./12)*AliHLTTPCTransform::GetPadPitchWidth(patch)*AliHLTTPCTransform::GetPadPitchWidth(patch);
1157 if(sy2 != 0){
1158 fpad2*=0.108; //constants are from offline studies
1159 if(patch<2)
1160 fpad2*=2.07;
1161 }
5e75f4e0 1162 }
6b15c309 1163 }
a912b63b 1164 // Float_t sz2=list[j].fTime2*list[j].fTotalCharge - list[j].fTime*list[j].fTime;
1165 Float_t sz2=(Float_t)list[j].fTime2*list[j].fTotalCharge - (Float_t)list[j].fTime*list[j].fTime;
1166 sz2/=q2;
1167 if(sz2 < 0){
1168 LOG(AliHLTTPCLog::kError,"AliHLTTPCClusterFinder::WriteClusters","Cluster width")
1169 <<"SigmaZ2 negative "<<sz2<<" on row "<<fCurrentRow<<" "<<fpad<<" "<<ftime<<ENDLOG;
1170 continue;
1171 } else {
5e75f4e0 1172 {
a912b63b 1173 ftime2 = (sz2 + 1./12)*AliHLTTPCTransform::GetZWidth()*AliHLTTPCTransform::GetZWidth();
1174 if(sz2 != 0) {
1175 ftime2 *= 0.169; //constants are from offline studies
1176 if(patch<2)
1177 ftime2 *= 1.77;
1178 }
5e75f4e0 1179 }
6b15c309 1180 }
1181 }
a912b63b 1182 if(fStdout==kTRUE)
1183 HLTInfo("WriteCluster: padrow %d pad %d +- %d time +- %d charge %d",fCurrentRow, fpad, fpad2, ftime, ftime2, list[j].fTotalCharge);
6b15c309 1184
5e75f4e0 1185 if (fFillRawClusters && fRawClusters.size()>(unsigned)counter) {
1186 fRawClusters[counter].SetPadRow(fCurrentRow);
1187 fRawClusters[counter].SetPad(fpad);
1188 fRawClusters[counter].SetTime(ftime);
1189 }
1190 {
a912b63b 1191 AliHLTTPCTransform::Slice2Sector(fCurrentSlice,fCurrentRow,thissector,thisrow);
5ee8e602 1192
1193 if(fOfflineTransform == NULL){
1194 AliHLTTPCTransform::Raw2Local(xyz,thissector,thisrow,fpad,ftime);
1195
1196 if(xyz[0]==0) LOG(AliHLTTPCLog::kError,"AliHLTTPCClustFinder","Cluster Finder")
1197 <<AliHLTTPCLog::kDec<<"Zero cluster"<<ENDLOG;
1198 if(fNClusters >= fMaxNClusters)
1199 {
1200 LOG(AliHLTTPCLog::kError,"AliHLTTPCClustFinder::WriteClusters","Cluster Finder")
1201 <<AliHLTTPCLog::kDec<<"Too many clusters "<<fNClusters<<ENDLOG;
1202 return;
1203 }
a912b63b 1204
5ee8e602 1205 fSpacePointData[counter].fX = xyz[0];
1206 // fSpacePointData[counter].fY = xyz[1];
1207 if(fCurrentSlice<18){
1208 fSpacePointData[counter].fY = xyz[1];
1209 }
1210 else{
1211 fSpacePointData[counter].fY = -1*xyz[1];
1212 }
1213 fSpacePointData[counter].fZ = xyz[2];
6b15c309 1214 }
a912b63b 1215 else{
2cb809a5 1216 Double_t x[3]={static_cast<Double_t>(thisrow),fpad+.5,ftime};
5ee8e602 1217 Int_t iSector[1]={thissector};
1218 fOfflineTransform->Transform(x,iSector,0,1);
a442ba54 1219 double y[3] = {x[0], x[1], x[2] };
1220
1221 if( fOfflineTPCParam && thissector<fOfflineTPCParam->GetNSector() ){
1222 TGeoHMatrix *alignment = fOfflineTPCParam->GetClusterMatrix( thissector );
1223 if ( alignment ) alignment->LocalToMaster( x, y);
1224 }
1225
1226 fSpacePointData[counter].fX = y[0];
1227 fSpacePointData[counter].fY = y[1];
1228 fSpacePointData[counter].fZ = y[2];
a912b63b 1229 }
8252a538 1230
5ee8e602 1231 }
a912b63b 1232
1233 fSpacePointData[counter].fCharge = list[j].fTotalCharge;
1234 fSpacePointData[counter].fPadRow = fCurrentRow;
1235 fSpacePointData[counter].fSigmaY2 = fpad2;
1236 fSpacePointData[counter].fSigmaZ2 = ftime2;
64defa03 1237
a912b63b 1238 fSpacePointData[counter].fQMax = list[j].fQMax;
2fdb1ae7 1239
1299d409 1240 fSpacePointData[counter].SetUsed(kFALSE); // only used / set in AliHLTTPCDisplay
1241 fSpacePointData[counter].SetTrackNumber(-1); // only used / set in AliHLTTPCDisplay
8252a538 1242
a912b63b 1243 Int_t patch=fCurrentPatch;
1244 if(patch==-1) patch=0; //never store negative patch number
a371a266 1245 fSpacePointData[counter].SetID( fCurrentSlice, patch, counter );
6b15c309 1246
5e75f4e0 1247 if (fFillRawClusters && fRawClusters.size()>(unsigned)counter) {
1248 fRawClusters[counter].SetSigmaY2(fSpacePointData[counter].fSigmaY2);
1249 fRawClusters[counter].SetSigmaZ2(fSpacePointData[counter].fSigmaZ2);
1250 fRawClusters[counter].SetCharge(fSpacePointData[counter].fCharge);
1251 fRawClusters[counter].SetQMax(fSpacePointData[counter].fQMax);
1252 }
1253
a912b63b 1254#ifdef do_mc
1255 Int_t trackID[3];
1256 GetTrackID((Int_t)rint(fpad),(Int_t)rint(ftime),trackID);
1257
1258 fSpacePointData[counter].fTrackID[0] = trackID[0];
1259 fSpacePointData[counter].fTrackID[1] = trackID[1];
1260 fSpacePointData[counter].fTrackID[2] = trackID[2];
1261
1262#endif
1263
1264 fNClusters++;
1265 counter++;
8252a538 1266 }
a912b63b 1267}
8252a538 1268
a912b63b 1269// STILL TO FIX ----------------------------------------------------------------------------
8252a538 1270
a912b63b 1271#ifdef do_mc
b783c3e8 1272void AliHLTTPCClusterFinder::GetTrackID(Int_t pad,Int_t time,Int_t *trackID) const {
a912b63b 1273 // see header file for class documentation
1274
1275 //get mc id
1276 AliHLTTPCDigitRowData *rowPt = (AliHLTTPCDigitRowData*)fDigitRowData;
8252a538 1277
a912b63b 1278 trackID[0]=trackID[1]=trackID[2]=-2;
1279 for(Int_t i=fFirstRow; i<=fLastRow; i++){
1280 if(rowPt->fRow < (UInt_t)fCurrentRow){
1281 AliHLTTPCMemHandler::UpdateRowPointer(rowPt);
1282 continue;
1283 }
1284 AliHLTTPCDigitData *digPt = (AliHLTTPCDigitData*)rowPt->fDigitData;
1285 for(UInt_t j=0; j<rowPt->fNDigit; j++){
1286 Int_t cpad = digPt[j].fPad;
1287 Int_t ctime = digPt[j].fTime;
1288 if(cpad != pad) continue;
1289 if(ctime != time) continue;
fec96a66 1290
a912b63b 1291 trackID[0] = digPt[j].fTrackID[0];
1292 trackID[1] = digPt[j].fTrackID[1];
1293 trackID[2] = digPt[j].fTrackID[2];
1294
1295 break;
1296 }
1297 break;
8252a538 1298 }
01f43166 1299}
a912b63b 1300#endif
64defa03 1301
a912b63b 1302
1303
1304void AliHLTTPCClusterFinder::WriteClusters(Int_t nclusters,AliHLTTPCClusters *list){//This is used when using the AliHLTTPCClusters class for cluster data
8252a538 1305 // see header file for class documentation
8252a538 1306
01f43166 1307 //write cluster to output pointer
1308 Int_t thisrow,thissector;
1309 UInt_t counter = fNClusters;
5e75f4e0 1310
1311 if (fFillRawClusters) {
1312 fRawClusters.resize(nclusters);
1313 }
01f43166 1314
1315 for(int j=0; j<nclusters; j++)
1316 {
1317 if(!list[j].fFlags) continue; //discard single pad clusters
1318 if(list[j].fTotalCharge < fThreshold) continue; //noise cluster
1319
1320 Float_t xyz[3];
1321 Float_t fpad =(Float_t)list[j].fPad / list[j].fTotalCharge;
1322 Float_t fpad2=fXYErr*fXYErr; //fixed given error
1323 Float_t ftime =(Float_t)list[j].fTime / list[j].fTotalCharge;
1324 Float_t ftime2=fZErr*fZErr; //fixed given error
1325
d8b2f74a 1326
01f43166 1327 if(fCalcerr) { //calc the errors, otherwice take the fixed error
1328 Int_t patch = AliHLTTPCTransform::GetPatch(fCurrentRow);
1329 UInt_t q2=list[j].fTotalCharge*list[j].fTotalCharge;
1330 Float_t sy2=list[j].fPad2 * list[j].fTotalCharge - list[j].fPad * list[j].fPad;
1331 sy2/=q2;
1332 if(sy2 < 0) {
1333 LOG(AliHLTTPCLog::kError,"AliHLTTPCClusterFinder::WriteClusters","Cluster width")
1334 <<"SigmaY2 negative "<<sy2<<" on row "<<fCurrentRow<<" "<<fpad<<" "<<ftime<<ENDLOG;
1335 continue;
1336 } else {
5e75f4e0 1337 {
01f43166 1338 fpad2 = (sy2 + 1./12)*AliHLTTPCTransform::GetPadPitchWidth(patch)*AliHLTTPCTransform::GetPadPitchWidth(patch);
1339 if(sy2 != 0){
1340 fpad2*=0.108; //constants are from offline studies
1341 if(patch<2)
1342 fpad2*=2.07;
1343 }
5e75f4e0 1344 }
01f43166 1345 }
1346 Float_t sz2=list[j].fTime2*list[j].fTotalCharge - list[j].fTime*list[j].fTime;
1347 sz2/=q2;
1348 if(sz2 < 0){
1349 LOG(AliHLTTPCLog::kError,"AliHLTTPCClusterFinder::WriteClusters","Cluster width")
1350 <<"SigmaZ2 negative "<<sz2<<" on row "<<fCurrentRow<<" "<<fpad<<" "<<ftime<<ENDLOG;
1351 continue;
1352 } else {
5e75f4e0 1353 {
01f43166 1354 ftime2 = (sz2 + 1./12)*AliHLTTPCTransform::GetZWidth()*AliHLTTPCTransform::GetZWidth();
1355 if(sz2 != 0) {
1356 ftime2 *= 0.169; //constants are from offline studies
1357 if(patch<2)
1358 ftime2 *= 1.77;
1359 }
5e75f4e0 1360 }
01f43166 1361 }
1362 }
1363 if(fStdout==kTRUE)
8252a538 1364 HLTInfo("WriteCluster: padrow %d pad %d +- %d time +- %d charge %d",fCurrentRow, fpad, fpad2, ftime, ftime2, list[j].fTotalCharge);
1365
5e75f4e0 1366 if (fFillRawClusters && fRawClusters.size()>(unsigned)counter) {
1367 fRawClusters[counter].SetPadRow(fCurrentRow);
1368 fRawClusters[counter].SetPad(fpad);
1369 fRawClusters[counter].SetTime(ftime);
1370 }
1371 {
01f43166 1372 AliHLTTPCTransform::Slice2Sector(fCurrentSlice,fCurrentRow,thissector,thisrow);
1373 AliHLTTPCTransform::Raw2Local(xyz,thissector,thisrow,fpad,ftime);
1374
1375 if(xyz[0]==0) LOG(AliHLTTPCLog::kError,"AliHLTTPCClustFinder","Cluster Finder")
1376 <<AliHLTTPCLog::kDec<<"Zero cluster"<<ENDLOG;
1377 if(fNClusters >= fMaxNClusters)
1378 {
1379 LOG(AliHLTTPCLog::kError,"AliHLTTPCClustFinder::WriteClusters","Cluster Finder")
1380 <<AliHLTTPCLog::kDec<<"Too many clusters "<<fNClusters<<ENDLOG;
1381 return;
1382 }
1383
1384 fSpacePointData[counter].fX = xyz[0];
6b15c309 1385 // fSpacePointData[counter].fY = xyz[1];
1386 if(fCurrentSlice<18){
1387 fSpacePointData[counter].fY = xyz[1];
1388 }
1389 else{
1390 fSpacePointData[counter].fY = -1*xyz[1];
1391 }
01f43166 1392 fSpacePointData[counter].fZ = xyz[2];
1393
01f43166 1394 }
1395
1396 fSpacePointData[counter].fCharge = list[j].fTotalCharge;
1397 fSpacePointData[counter].fPadRow = fCurrentRow;
1398 fSpacePointData[counter].fSigmaY2 = fpad2;
1399 fSpacePointData[counter].fSigmaZ2 = ftime2;
1400
a912b63b 1401 fSpacePointData[counter].fQMax = list[j].fQMax;
6b15c309 1402
1299d409 1403 fSpacePointData[counter].SetUsed(kFALSE); // only used / set in AliHLTTPCDisplay
1404 fSpacePointData[counter].SetTrackNumber(-1); // only used / set in AliHLTTPCDisplay
01f43166 1405
1406 Int_t patch=fCurrentPatch;
1407 if(patch==-1) patch=0; //never store negative patch number
a371a266 1408 fSpacePointData[counter].SetID( fCurrentSlice, patch, counter );
01f43166 1409
5e75f4e0 1410 if (fFillRawClusters && fRawClusters.size()>(unsigned)counter) {
1411 fRawClusters[counter].SetSigmaY2(fSpacePointData[counter].fSigmaY2);
1412 fRawClusters[counter].SetSigmaZ2(fSpacePointData[counter].fSigmaZ2);
1413 fRawClusters[counter].SetCharge(fSpacePointData[counter].fCharge);
1414 fRawClusters[counter].SetQMax(fSpacePointData[counter].fQMax);
1415 }
1416
01f43166 1417#ifdef do_mc
1418 Int_t trackID[3];
1419 GetTrackID((Int_t)rint(fpad),(Int_t)rint(ftime),trackID);
1420
1421 fSpacePointData[counter].fTrackID[0] = trackID[0];
1422 fSpacePointData[counter].fTrackID[1] = trackID[1];
1423 fSpacePointData[counter].fTrackID[2] = trackID[2];
1424
01f43166 1425#endif
1426
1427 fNClusters++;
1428 counter++;
1429 }
1430}