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