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