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