]> git.uio.no Git - u/mrichter/AliRoot.git/blame - HLT/TPCLib/AliHLTTPCClusterFinder.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Bug fix for run ranges having gaps in the actual data.
[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){
1022 fClustersMCInfo.erase(fClustersMCInfo.begin()+j); // remove the mc ifo for this cluster since it is not taken into account
1023 }
1024 continue; //discard single pad clusters
1025 }
1026 if(list[j].fTotalCharge < fThreshold){
1027 if(fDoMC){
1028 fClustersMCInfo.erase(fClustersMCInfo.begin()+j); // remove the mc ifo for this cluster since it is not taken into account
1029 }
1030 continue; //noise cluster
1031 }
a38a7850 1032 Float_t xyz[3];
1033 Float_t fpad =(Float_t)list[j].fPad / list[j].fTotalCharge;
1034 Float_t fpad2=fXYErr*fXYErr; //fixed given error
1035 Float_t ftime =(Float_t)list[j].fTime / list[j].fTotalCharge;
1036 Float_t ftime2=fZErr*fZErr; //fixed given error
1037
738c049f 1038
aff6e981 1039
a912b63b 1040 if(fUnsorted){
1041 fCurrentRow=list[j].fRow;
8252a538 1042 }
8252a538 1043
a912b63b 1044
1045 if(fCalcerr) { //calc the errors, otherwice take the fixed error
1046 Int_t patch = AliHLTTPCTransform::GetPatch(fCurrentRow);
1047 UInt_t q2=list[j].fTotalCharge*list[j].fTotalCharge;
1048 // Float_t sy2=list[j].fPad2 * list[j].fTotalCharge - list[j].fPad * list[j].fPad;
1049 Float_t sy2=(Float_t)list[j].fPad2 * list[j].fTotalCharge - (Float_t)list[j].fPad * list[j].fPad;
1050 sy2/=q2;
1051 if(sy2 < 0) {
1052 LOG(AliHLTTPCLog::kError,"AliHLTTPCClusterFinder::WriteClusters","Cluster width")
1053 <<"SigmaY2 negative "<<sy2<<" on row "<<fCurrentRow<<" "<<fpad<<" "<<ftime<<ENDLOG;
1054 continue;
1055 } else {
1056 if(!fRawSP){
1057 fpad2 = (sy2 + 1./12)*AliHLTTPCTransform::GetPadPitchWidth(patch)*AliHLTTPCTransform::GetPadPitchWidth(patch);
1058 if(sy2 != 0){
1059 fpad2*=0.108; //constants are from offline studies
1060 if(patch<2)
1061 fpad2*=2.07;
1062 }
1063 } else fpad2=sy2; //take the width not the error
6b15c309 1064 }
a912b63b 1065 // Float_t sz2=list[j].fTime2*list[j].fTotalCharge - list[j].fTime*list[j].fTime;
1066 Float_t sz2=(Float_t)list[j].fTime2*list[j].fTotalCharge - (Float_t)list[j].fTime*list[j].fTime;
1067 sz2/=q2;
1068 if(sz2 < 0){
1069 LOG(AliHLTTPCLog::kError,"AliHLTTPCClusterFinder::WriteClusters","Cluster width")
1070 <<"SigmaZ2 negative "<<sz2<<" on row "<<fCurrentRow<<" "<<fpad<<" "<<ftime<<ENDLOG;
1071 continue;
1072 } else {
1073 if(!fRawSP){
1074 ftime2 = (sz2 + 1./12)*AliHLTTPCTransform::GetZWidth()*AliHLTTPCTransform::GetZWidth();
1075 if(sz2 != 0) {
1076 ftime2 *= 0.169; //constants are from offline studies
1077 if(patch<2)
1078 ftime2 *= 1.77;
1079 }
1080 } else ftime2=sz2; //take the width, not the error
6b15c309 1081 }
1082 }
a912b63b 1083 if(fStdout==kTRUE)
1084 HLTInfo("WriteCluster: padrow %d pad %d +- %d time +- %d charge %d",fCurrentRow, fpad, fpad2, ftime, ftime2, list[j].fTotalCharge);
6b15c309 1085
a912b63b 1086 if(!fRawSP){
1087 AliHLTTPCTransform::Slice2Sector(fCurrentSlice,fCurrentRow,thissector,thisrow);
1088 AliHLTTPCTransform::Raw2Local(xyz,thissector,thisrow,fpad,ftime);
1089
1090 if(xyz[0]==0) LOG(AliHLTTPCLog::kError,"AliHLTTPCClustFinder","Cluster Finder")
1091 <<AliHLTTPCLog::kDec<<"Zero cluster"<<ENDLOG;
1092 if(fNClusters >= fMaxNClusters)
1093 {
1094 LOG(AliHLTTPCLog::kError,"AliHLTTPCClustFinder::WriteClusters","Cluster Finder")
1095 <<AliHLTTPCLog::kDec<<"Too many clusters "<<fNClusters<<ENDLOG;
1096 return;
1097 }
1098
1099 fSpacePointData[counter].fX = xyz[0];
1100 // fSpacePointData[counter].fY = xyz[1];
1101 if(fCurrentSlice<18){
1102 fSpacePointData[counter].fY = xyz[1];
6b15c309 1103 }
a912b63b 1104 else{
1105 fSpacePointData[counter].fY = -1*xyz[1];
1106 }
1107 fSpacePointData[counter].fZ = xyz[2];
8252a538 1108
a912b63b 1109 } else {
1110 fSpacePointData[counter].fX = fCurrentRow;
1111 fSpacePointData[counter].fY = fpad;
1112 fSpacePointData[counter].fZ = ftime;
8252a538 1113 }
a912b63b 1114
1115 fSpacePointData[counter].fCharge = list[j].fTotalCharge;
1116 fSpacePointData[counter].fPadRow = fCurrentRow;
1117 fSpacePointData[counter].fSigmaY2 = fpad2;
1118 fSpacePointData[counter].fSigmaZ2 = ftime2;
64defa03 1119
a912b63b 1120 fSpacePointData[counter].fQMax = list[j].fQMax;
2fdb1ae7 1121
a912b63b 1122 fSpacePointData[counter].fUsed = kFALSE; // only used / set in AliHLTTPCDisplay
1123 fSpacePointData[counter].fTrackN = -1; // only used / set in AliHLTTPCDisplay
8252a538 1124
a912b63b 1125 Int_t patch=fCurrentPatch;
1126 if(patch==-1) patch=0; //never store negative patch number
1127 fSpacePointData[counter].fID = counter
1128 +((fCurrentSlice&0x7f)<<25)+((patch&0x7)<<22);//Uli
6b15c309 1129
a912b63b 1130#ifdef do_mc
1131 Int_t trackID[3];
1132 GetTrackID((Int_t)rint(fpad),(Int_t)rint(ftime),trackID);
1133
1134 fSpacePointData[counter].fTrackID[0] = trackID[0];
1135 fSpacePointData[counter].fTrackID[1] = trackID[1];
1136 fSpacePointData[counter].fTrackID[2] = trackID[2];
1137
1138#endif
1139
1140 fNClusters++;
1141 counter++;
8252a538 1142 }
a912b63b 1143}
8252a538 1144
a912b63b 1145// STILL TO FIX ----------------------------------------------------------------------------
8252a538 1146
a912b63b 1147#ifdef do_mc
1148void AliHLTTPCClusterFinder::GetTrackID(Int_t pad,Int_t time,Int_t *trackID){
1149 // see header file for class documentation
1150
1151 //get mc id
1152 AliHLTTPCDigitRowData *rowPt = (AliHLTTPCDigitRowData*)fDigitRowData;
8252a538 1153
a912b63b 1154 trackID[0]=trackID[1]=trackID[2]=-2;
1155 for(Int_t i=fFirstRow; i<=fLastRow; i++){
1156 if(rowPt->fRow < (UInt_t)fCurrentRow){
1157 AliHLTTPCMemHandler::UpdateRowPointer(rowPt);
1158 continue;
1159 }
1160 AliHLTTPCDigitData *digPt = (AliHLTTPCDigitData*)rowPt->fDigitData;
1161 for(UInt_t j=0; j<rowPt->fNDigit; j++){
1162 Int_t cpad = digPt[j].fPad;
1163 Int_t ctime = digPt[j].fTime;
1164 if(cpad != pad) continue;
1165 if(ctime != time) continue;
fec96a66 1166
a912b63b 1167 trackID[0] = digPt[j].fTrackID[0];
1168 trackID[1] = digPt[j].fTrackID[1];
1169 trackID[2] = digPt[j].fTrackID[2];
1170
1171 break;
1172 }
1173 break;
8252a538 1174 }
01f43166 1175}
a912b63b 1176#endif
64defa03 1177
a912b63b 1178
1179
1180void AliHLTTPCClusterFinder::WriteClusters(Int_t nclusters,AliHLTTPCClusters *list){//This is used when using the AliHLTTPCClusters class for cluster data
8252a538 1181 // see header file for class documentation
8252a538 1182
01f43166 1183 //write cluster to output pointer
1184 Int_t thisrow,thissector;
1185 UInt_t counter = fNClusters;
1186
1187 for(int j=0; j<nclusters; j++)
1188 {
1189 if(!list[j].fFlags) continue; //discard single pad clusters
1190 if(list[j].fTotalCharge < fThreshold) continue; //noise cluster
1191
1192 Float_t xyz[3];
1193 Float_t fpad =(Float_t)list[j].fPad / list[j].fTotalCharge;
1194 Float_t fpad2=fXYErr*fXYErr; //fixed given error
1195 Float_t ftime =(Float_t)list[j].fTime / list[j].fTotalCharge;
1196 Float_t ftime2=fZErr*fZErr; //fixed given error
1197
d8b2f74a 1198
01f43166 1199 if(fCalcerr) { //calc the errors, otherwice take the fixed error
1200 Int_t patch = AliHLTTPCTransform::GetPatch(fCurrentRow);
1201 UInt_t q2=list[j].fTotalCharge*list[j].fTotalCharge;
1202 Float_t sy2=list[j].fPad2 * list[j].fTotalCharge - list[j].fPad * list[j].fPad;
1203 sy2/=q2;
1204 if(sy2 < 0) {
1205 LOG(AliHLTTPCLog::kError,"AliHLTTPCClusterFinder::WriteClusters","Cluster width")
1206 <<"SigmaY2 negative "<<sy2<<" on row "<<fCurrentRow<<" "<<fpad<<" "<<ftime<<ENDLOG;
1207 continue;
1208 } else {
1209 if(!fRawSP){
1210 fpad2 = (sy2 + 1./12)*AliHLTTPCTransform::GetPadPitchWidth(patch)*AliHLTTPCTransform::GetPadPitchWidth(patch);
1211 if(sy2 != 0){
1212 fpad2*=0.108; //constants are from offline studies
1213 if(patch<2)
1214 fpad2*=2.07;
1215 }
1216 } else fpad2=sy2; //take the width not the error
1217 }
1218 Float_t sz2=list[j].fTime2*list[j].fTotalCharge - list[j].fTime*list[j].fTime;
1219 sz2/=q2;
1220 if(sz2 < 0){
1221 LOG(AliHLTTPCLog::kError,"AliHLTTPCClusterFinder::WriteClusters","Cluster width")
1222 <<"SigmaZ2 negative "<<sz2<<" on row "<<fCurrentRow<<" "<<fpad<<" "<<ftime<<ENDLOG;
1223 continue;
1224 } else {
1225 if(!fRawSP){
1226 ftime2 = (sz2 + 1./12)*AliHLTTPCTransform::GetZWidth()*AliHLTTPCTransform::GetZWidth();
1227 if(sz2 != 0) {
1228 ftime2 *= 0.169; //constants are from offline studies
1229 if(patch<2)
1230 ftime2 *= 1.77;
1231 }
1232 } else ftime2=sz2; //take the width, not the error
1233 }
1234 }
1235 if(fStdout==kTRUE)
8252a538 1236 HLTInfo("WriteCluster: padrow %d pad %d +- %d time +- %d charge %d",fCurrentRow, fpad, fpad2, ftime, ftime2, list[j].fTotalCharge);
1237
01f43166 1238 if(!fRawSP){
1239 AliHLTTPCTransform::Slice2Sector(fCurrentSlice,fCurrentRow,thissector,thisrow);
1240 AliHLTTPCTransform::Raw2Local(xyz,thissector,thisrow,fpad,ftime);
1241
1242 if(xyz[0]==0) LOG(AliHLTTPCLog::kError,"AliHLTTPCClustFinder","Cluster Finder")
1243 <<AliHLTTPCLog::kDec<<"Zero cluster"<<ENDLOG;
1244 if(fNClusters >= fMaxNClusters)
1245 {
1246 LOG(AliHLTTPCLog::kError,"AliHLTTPCClustFinder::WriteClusters","Cluster Finder")
1247 <<AliHLTTPCLog::kDec<<"Too many clusters "<<fNClusters<<ENDLOG;
1248 return;
1249 }
1250
1251 fSpacePointData[counter].fX = xyz[0];
6b15c309 1252 // fSpacePointData[counter].fY = xyz[1];
1253 if(fCurrentSlice<18){
1254 fSpacePointData[counter].fY = xyz[1];
1255 }
1256 else{
1257 fSpacePointData[counter].fY = -1*xyz[1];
1258 }
01f43166 1259 fSpacePointData[counter].fZ = xyz[2];
1260
1261 } else {
1262 fSpacePointData[counter].fX = fCurrentRow;
1263 fSpacePointData[counter].fY = fpad;
1264 fSpacePointData[counter].fZ = ftime;
1265 }
1266
1267 fSpacePointData[counter].fCharge = list[j].fTotalCharge;
1268 fSpacePointData[counter].fPadRow = fCurrentRow;
1269 fSpacePointData[counter].fSigmaY2 = fpad2;
1270 fSpacePointData[counter].fSigmaZ2 = ftime2;
1271
a912b63b 1272 fSpacePointData[counter].fQMax = list[j].fQMax;
6b15c309 1273
01f43166 1274 fSpacePointData[counter].fUsed = kFALSE; // only used / set in AliHLTTPCDisplay
1275 fSpacePointData[counter].fTrackN = -1; // only used / set in AliHLTTPCDisplay
1276
1277 Int_t patch=fCurrentPatch;
1278 if(patch==-1) patch=0; //never store negative patch number
1279 fSpacePointData[counter].fID = counter
1280 +((fCurrentSlice&0x7f)<<25)+((patch&0x7)<<22);//Uli
1281
1282#ifdef do_mc
1283 Int_t trackID[3];
1284 GetTrackID((Int_t)rint(fpad),(Int_t)rint(ftime),trackID);
1285
1286 fSpacePointData[counter].fTrackID[0] = trackID[0];
1287 fSpacePointData[counter].fTrackID[1] = trackID[1];
1288 fSpacePointData[counter].fTrackID[2] = trackID[2];
1289
01f43166 1290#endif
1291
1292 fNClusters++;
1293 counter++;
1294 }
1295}