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