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