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