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