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