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